Organize data

Merge output data from multiple batches, select metabolite columns, add participant info and other data

acg_batch1 <- read.csv("/Users/meaghan/OneDrive - University of Calgary/Preschool_data/MRS/data/ACG_LCM_CorrectedTable_12-14-2021-16-15.csv")

acg_batch2 <- read.csv("/Users/meaghan/OneDrive - University of Calgary/Preschool_data/MRS/data/ACG_LCM_CorrectedTable_12-14-2021-16-21.csv")

acg_LCM_corr <- bind_rows(acg_batch1, acg_batch2)
anyDuplicated(acg_LCM_corr$File_ID) #check for duplicated values
## [1] 106
acg_LCM_corr<-unique(acg_LCM_corr) #remove duplicate rows
acg_LCM_corr["4", "Cho_GPC_PCh_conc_molar"]<-NA #Cho conc invalid for PS10008_PS14_046_2386 due to spike in spectrum
acg_LCM_corr["4", "Cho_GPC_PCh_conc_molal"]<-NA #Cho conc invalid for PS10008_PS14_046_2386 due to spike in spectrum

write.csv(acg_LCM_corr, "/Users/meaghan/OneDrive - University of Calgary/Preschool_data/MRS/data/ACG_LCM_CorrectedTable_12-14-2021_all.csv", row.names = FALSE)

acg_molal <- acg_LCM_corr %>%
  dplyr::select(File_ID, fGM, fWM, fCSF, NAA_conc_molal, Cre_PCr_conc_molal, Cho_GPC_PCh_conc_molal, Ins_conc_molal, Glu_conc_molal, Glu_Gln_conc_molal)

acg_demo <- read.csv("/Users/meaghan/OneDrive - University of Calgary/Preschool_data/MRS/data/acg_analysis_subs.csv")

acg_dat <- left_join(acg_demo, acg_molal, by = "File_ID")
write.csv(acg_dat, "/Users/meaghan/OneDrive - University of Calgary/Preschool_data/MRS/data/acg_data_molal_n113.csv", row.names = FALSE)

lag_batch1 <- read.csv("/Users/meaghan/OneDrive - University of Calgary/Preschool_data/MRS/data/LAG_LCM_CorrectedTable_12-13-2021-15-18.csv")

lag_batch2 <- read.csv("/Users/meaghan/OneDrive - University of Calgary/Preschool_data/MRS/data/LAG_LCM_CorrectedTable_12-14-2021-11-41.csv")

lag_batch3 <- read.csv("/Users/meaghan/OneDrive - University of Calgary/Preschool_data/MRS/data/LAG_LCM_CorrectedTable_12-14-2021-12-18.csv")

lag_LCM_corr <- bind_rows(lag_batch1, lag_batch2, lag_batch3)
write.csv(lag_LCM_corr, "/Users/meaghan/OneDrive - University of Calgary/Preschool_data/MRS/data/LAG_LCM_CorrectedTable_12-14-2021_all.csv", row.names = FALSE)

lag_molal <- lag_LCM_corr %>%
  dplyr::select(File_ID, fGM, fWM, fCSF, NAA_conc_molal, Cre_PCr_conc_molal, Cho_GPC_PCh_conc_molal, Ins_conc_molal, Glu_conc_molal, Glu_Gln_conc_molal)

lag_demo <- read.csv("/Users/meaghan/OneDrive - University of Calgary/Preschool_data/MRS/data/lag_analysis_subs.csv")

lag_dat <- left_join(lag_demo, lag_molal, by = "File_ID")
lag_dat$subj_id<-as.character(lag_dat$subj_id) #make subj_id character type

write.csv(lag_dat, "/Users/meaghan/OneDrive - University of Calgary/Preschool_data/MRS/data/lag_data_molal_n321.csv", row.names = FALSE)

##Check data distribution & outliers Create plots to check data distribution and identify outliers

ACG

plot_age <- ggplot(acg_dat, aes(mri_age_y))
plot_age + geom_density()

#age skewed

plot_naa <- ggplot(acg_dat, aes(NAA_conc_molal))
plot_naa + geom_density()

plot_naa + geom_boxplot()

rosnerTest(acg_dat$NAA_conc_molal, k = 4)
## $distribution
## [1] "Normal"
## 
## $statistic
##      R.1      R.2      R.3      R.4 
## 2.759147 2.825357 2.694361 2.721635 
## 
## $sample.size
## [1] 113
## 
## $parameters
## k 
## 4 
## 
## $alpha
## [1] 0.05
## 
## $crit.value
## lambda.1 lambda.2 lambda.3 lambda.4 
## 3.425263 3.422302 3.419309 3.416284 
## 
## $n.outliers
## [1] 0
## 
## $alternative
## [1] "Up to 4 observations are not\n                                 from the same Distribution."
## 
## $method
## [1] "Rosner's Test for Outliers"
## 
## $data
##   [1] 13.36240 15.05227 14.57780 13.83339 14.71206 16.20111 15.44740 16.98446
##   [9] 15.67072 15.95324 14.87793 12.74216 16.29728 14.71133 14.66048 13.28760
##  [17] 14.59532 13.55773 16.85293 15.33699 12.45017 15.04043 15.27049 11.24177
##  [25] 14.98401 14.18774 15.83415 13.91788 17.10858 12.09629 14.70242 14.40906
##  [33] 16.46456 15.58993 15.03289 14.84971 11.30097 13.87568 15.24533 15.99486
##  [41] 14.53574 12.16795 14.20751 14.18561 16.47436 14.52629 15.00648 16.20561
##  [49] 13.22605 13.57437 14.05939 13.81423 17.33147 16.31303 16.29984 16.85100
##  [57] 15.24641 14.93994 15.14031 15.13625 18.29604 13.89531 13.42919 15.16179
##  [65] 13.87656 15.48760 12.64660 16.01636 13.64857 16.43270 15.59565 11.64927
##  [73] 15.46789 14.18371 12.72496 16.73666 16.33378 15.15593 15.86369 15.28939
##  [81] 15.31529 14.67805 13.38081 15.01442 13.01089 17.11696 14.54570 15.04389
##  [89] 14.64002 15.57552 15.21758 15.80866 15.16132 14.92979 14.18939 16.96225
##  [97] 13.86479 15.05179 15.22074 15.65731 15.51013 14.39392 15.03669 16.39137
## [105] 14.47850 14.65614 15.94476 16.78988 16.07591 13.54665 13.48426 14.73463
## [113] 15.68117
## 
## $data.name
## [1] "acg_dat$NAA_conc_molal"
## 
## $bad.obs
## [1] 0
## 
## $all.stats
##   i   Mean.i     SD.i    Value Obs.Num    R.i+1 lambda.i+1 Outlier
## 1 0 14.88986 1.322180 11.24177      24 2.759147   3.425263   FALSE
## 2 1 14.92243 1.281772 11.30097      37 2.825357   3.422302   FALSE
## 3 2 14.95506 1.239990 18.29604      61 2.694361   3.419309   FALSE
## 4 3 14.92469 1.203472 11.64927      72 2.721635   3.416284   FALSE
## 
## attr(,"class")
## [1] "gofOutlier"
plot_cre <- ggplot(acg_dat, aes(Cre_PCr_conc_molal))
plot_cre + geom_density()

plot_cre + geom_boxplot()

rosnerTest(acg_dat$Cre_PCr_conc_molal, k = 5)
## $distribution
## [1] "Normal"
## 
## $statistic
##      R.1      R.2      R.3      R.4      R.5 
## 3.394290 2.833607 2.707716 2.578893 2.649738 
## 
## $sample.size
## [1] 113
## 
## $parameters
## k 
## 5 
## 
## $alpha
## [1] 0.05
## 
## $crit.value
## lambda.1 lambda.2 lambda.3 lambda.4 lambda.5 
## 3.425263 3.422302 3.419309 3.416284 3.413225 
## 
## $n.outliers
## [1] 0
## 
## $alternative
## [1] "Up to 5 observations are not\n                                 from the same Distribution."
## 
## $method
## [1] "Rosner's Test for Outliers"
## 
## $data
##   [1] 10.869437 10.445542 14.138594 11.979460 11.842544 12.567296 12.183610
##   [8] 13.315617 13.532807 13.910100 12.102060  9.609688 12.566391 11.689508
##  [15] 12.650725 10.976876 11.888348 11.823341 12.975541 12.087174 10.745566
##  [22] 12.005758 13.958950 11.012768 13.437536 12.359373 12.418194 10.775703
##  [29] 13.584172 11.251205 11.993304 12.190118 14.756858 13.780875 12.273252
##  [36] 11.243397  8.406874 11.279114 11.602557 13.761404 12.044563 10.276615
##  [43] 10.918920 10.896575 11.867174 12.189638 12.580644 12.364280 11.517322
##  [50] 10.442494 11.903800 12.566843 12.411654 12.733123 12.852863 14.066514
##  [57] 10.742053 11.977272 12.534036 12.467809 13.778825 11.508522 11.740032
##  [64] 11.569008 10.891687 12.474281 11.395342 12.716648 12.033736 13.646045
##  [71] 12.873670  9.632426 11.901341 12.610355  9.187364 13.087309 12.160354
##  [78] 11.717870 12.392428 11.411606 12.417766 11.214681 11.894967 11.689931
##  [85] 11.069854 13.648131 12.837410 11.764882 12.329981 11.666319 12.825568
##  [92] 13.032058 12.630243 12.355787 11.038718  9.731136 12.044142 11.736949
##  [99] 11.772225 11.427900 12.588592 11.074387 11.520514 12.986814 12.563423
## [106] 12.334515 13.095183 12.401466 11.143206 11.414014 11.499070 12.331493
## [113] 11.900989
## 
## $data.name
## [1] "acg_dat$Cre_PCr_conc_molal"
## 
## $bad.obs
## [1] 0
## 
## $all.stats
##   i   Mean.i      SD.i     Value Obs.Num    R.i+1 lambda.i+1 Outlier
## 1 0 12.03590 1.0691568  8.406874      37 3.394290   3.425263   FALSE
## 2 1 12.06830 1.0167044  9.187364      75 2.833607   3.422302   FALSE
## 3 2 12.09426 0.9833373 14.756858      33 2.707716   3.419309   FALSE
## 4 3 12.07005 0.9540396  9.609688      12 2.578893   3.416284   FALSE
## 5 4 12.09263 0.9284692  9.632426      72 2.649738   3.413225   FALSE
## 
## attr(,"class")
## [1] "gofOutlier"
plot_cho <- ggplot(acg_dat, aes(Cho_GPC_PCh_conc_molal))
plot_cho + geom_density()
## Warning: Removed 1 rows containing non-finite values (stat_density).

plot_cho + geom_boxplot()
## Warning: Removed 1 rows containing non-finite values (stat_boxplot).

rosnerTest(acg_dat$Cho_GPC_PCh_conc_molal, k = 4)
## Warning in rosnerTest(acg_dat$Cho_GPC_PCh_conc_molal, k = 4): 1 observations
## with NA/NaN/Inf in 'x' removed.
## $distribution
## [1] "Normal"
## 
## $statistic
##      R.1      R.2      R.3      R.4 
## 2.692059 2.486294 2.542571 2.587426 
## 
## $sample.size
## [1] 112
## 
## $parameters
## k 
## 4 
## 
## $alpha
## [1] 0.05
## 
## $crit.value
## lambda.1 lambda.2 lambda.3 lambda.4 
## 3.422302 3.419309 3.416284 3.413225 
## 
## $n.outliers
## [1] 0
## 
## $alternative
## [1] "Up to 4 observations are not\n                                 from the same Distribution."
## 
## $method
## [1] "Rosner's Test for Outliers"
## 
## $data
##   [1] 3.054728 2.338141 2.483788 2.615586 2.904744 3.029351 2.945605 2.946224
##   [9] 2.469310 2.647417 2.329234 2.525650 2.964839 2.377024 2.316960 2.388015
##  [17] 2.587590 2.964150 2.987201 1.884924 2.294818 2.160623 1.723047 2.273437
##  [25] 2.486538 2.364455 2.110030 3.042283 2.145223 2.626855 2.657139 3.307897
##  [33] 2.993478 2.528009 2.694001 2.751861 2.169544 2.576529 2.777803 2.298135
##  [41] 2.093485 2.245856 1.887058 3.332061 2.683945 2.960541 2.306091 2.400287
##  [49] 2.262987 2.373107 2.418340 2.428692 2.651362 2.724707 2.627568 3.213338
##  [57] 2.541659 2.436647 2.503022 3.168533 2.273151 2.230455 2.575725 2.463493
##  [65] 2.964749 2.427496 2.382196 2.398624 2.592440 2.075228 1.967159 2.210616
##  [73] 2.178031 1.572637 2.265866 2.193475 2.435427 2.749229 2.738993 2.096637
##  [81] 2.376787 2.150157 2.031053 1.736687 2.549619 2.032952 2.153716 2.354125
##  [89] 2.585598 2.624180 2.597158 2.503254 2.679817 2.109408 2.755517 2.618999
##  [97] 2.254582 2.379855 2.105926 2.389925 2.621861 2.349337 2.962557 2.620008
## [105] 2.643999 2.559238 3.322683 2.403398 2.301871 2.505683 2.762249 2.704160
## 
## $data.name
## [1] "acg_dat$Cho_GPC_PCh_conc_molal"
## 
## $bad.obs
## [1] 1
## 
## $all.stats
##   i   Mean.i      SD.i    Value Obs.Num    R.i+1 lambda.i+1 Outlier
## 1 0 2.495888 0.3429536 1.572637      75 2.692059   3.422302   FALSE
## 2 1 2.504206 0.3329678 3.332061      45 2.486294   3.419309   FALSE
## 3 2 2.496680 0.3248691 3.322683     108 2.542571   3.416284   FALSE
## 4 3 2.489102 0.3164519 3.307897      33 2.587426   3.413225   FALSE
## 
## attr(,"class")
## [1] "gofOutlier"
plot_ins <- ggplot(acg_dat, aes(Ins_conc_molal))
plot_ins + geom_density()

shapiro.test(acg_dat$Ins_conc_molal) #not normal
## 
##  Shapiro-Wilk normality test
## 
## data:  acg_dat$Ins_conc_molal
## W = 0.96211, p-value = 0.002754
plot_ins + geom_boxplot()

rosnerTest(acg_dat$Ins_conc_molal, k = 3) #outlier value 1.471244, observation 84 (PS18_013) - spectrum showed issue, exclude from analysis 
## $distribution
## [1] "Normal"
## 
## $statistic
##      R.1      R.2      R.3 
## 3.753416 3.044570 3.051178 
## 
## $sample.size
## [1] 113
## 
## $parameters
## k 
## 3 
## 
## $alpha
## [1] 0.05
## 
## $crit.value
## lambda.1 lambda.2 lambda.3 
## 3.425263 3.422302 3.419309 
## 
## $n.outliers
## [1] 1
## 
## $alternative
## [1] "Up to 3 observations are not\n                                 from the same Distribution."
## 
## $method
## [1] "Rosner's Test for Outliers"
## 
## $data
##   [1]  7.633881  8.178435  7.218240  7.696923  8.696692  8.600324  9.402207
##   [8]  8.875539  6.779586  5.935284 10.304907  5.028532  7.058793  9.067489
##  [15]  6.700483  7.361395  8.124760  8.682546  6.553380  6.320046  7.315123
##  [22]  7.446418  7.281237  7.102747  9.505161  7.839894  7.887706  7.616099
##  [29]  7.601724  6.347168  7.912922  5.656884  9.656654  9.672370  8.341248
##  [36]  8.816154  5.196891  6.304340  6.789199  8.018694  7.594973  7.416233
##  [43]  7.406756  3.198494  2.985922  9.898020  8.079535  8.703354  7.157223
##  [50]  6.213808  8.399767  9.174211  9.656691  9.335033 10.080614  4.813291
##  [57]  9.728369  7.521052  7.441993  8.051283 10.887902  8.695335  5.872553
##  [64]  6.418561  5.242212  8.749670  7.765345  8.185832  7.950928  8.394652
##  [71]  9.583041  6.202420  7.829477  6.209769  4.285017  8.677590  7.775045
##  [78]  9.334317  8.684385  7.236647  7.742767  6.219638  9.105159  1.471244
##  [85]  5.218826  9.407651  7.634559  6.064850  8.785683  3.770173  7.618745
##  [92]  8.763466  7.123140  6.399259  6.804541  5.855502  9.364210  6.203701
##  [99] 10.565682  9.856594  7.453487  7.728481  7.846516  9.947994  8.640216
## [106]  7.758280  8.476271  8.672477  6.437713  5.739577  6.673856  7.596796
## [113]  7.229168
## 
## $data.name
## [1] "acg_dat$Ins_conc_molal"
## 
## $bad.obs
## [1] 0
## 
## $all.stats
##   i   Mean.i     SD.i    Value Obs.Num    R.i+1 lambda.i+1 Outlier
## 1 0 7.588864 1.629881 1.471244      84 3.753416   3.425263    TRUE
## 2 1 7.643485 1.529793 2.985922      45 3.044570   3.422302   FALSE
## 3 2 7.685445 1.470564 3.198494      44 3.051178   3.419309   FALSE
## 
## attr(,"class")
## [1] "gofOutlier"
acg_dat["84", "Ins_conc_molal"] <- NA #replace outlier value with NA to exclude from analysis

plot_glx <- ggplot(acg_dat, aes(Glu_Gln_conc_molal))
plot_glx + geom_density()

plot_glx + geom_boxplot()

rosnerTest(acg_dat$Glu_Gln_conc_molal, k = 1)
## $distribution
## [1] "Normal"
## 
## $statistic
##      R.1 
## 2.840029 
## 
## $sample.size
## [1] 113
## 
## $parameters
## k 
## 1 
## 
## $alpha
## [1] 0.05
## 
## $crit.value
## lambda.1 
## 3.425263 
## 
## $n.outliers
## [1] 0
## 
## $alternative
## [1] "Up to 1 observations are not\n                                 from the same Distribution."
## 
## $method
## [1] "Rosner's Test for Outliers"
## 
## $data
##   [1] 22.16208 24.54895 30.79902 25.15445 23.12144 23.45996 27.95219 28.85494
##   [9] 28.47293 30.83024 24.02170 20.83895 26.81225 25.81936 25.83210 28.15351
##  [17] 28.69984 24.45281 25.87441 26.10568 24.90431 27.69072 29.16664 25.73952
##  [25] 31.00567 25.29233 29.35298 28.81394 30.45098 26.42377 25.25458 30.67838
##  [33] 29.65833 34.25604 26.24353 28.06087 21.93668 25.00839 27.05607 29.92424
##  [41] 27.12366 23.65348 26.82576 24.86136 26.92752 23.98651 26.21693 28.06806
##  [49] 26.26637 27.44238 23.98526 25.40089 22.59948 23.46421 28.23542 33.70878
##  [57] 26.99835 25.59893 26.95525 25.77541 32.34599 24.66776 25.46001 25.71687
##  [65] 29.73266 25.54233 27.95401 29.18637 26.09421 27.74434 27.70571 21.88957
##  [73] 26.72328 30.44184 21.32831 31.44932 28.73199 25.08766 29.57015 24.10930
##  [81] 27.84161 24.85586 29.21009 20.93551 23.34771 27.38610 26.35670 26.10330
##  [89] 28.90132 26.67329 24.94565 29.44010 29.32746 26.23233 24.32566 24.39736
##  [97] 26.88103 30.91320 28.56651 27.47953 25.89977 24.22166 25.89081 28.12079
## [105] 29.77257 28.04112 29.21096 29.61668 27.08694 28.10033 24.05989 26.66364
## [113] 26.45077
## 
## $data.name
## [1] "acg_dat$Glu_Gln_conc_molal"
## 
## $bad.obs
## [1] 0
## 
## $all.stats
##   i   Mean.i     SD.i    Value Obs.Num    R.i+1 lambda.i+1 Outlier
## 1 0 26.82935 2.615006 34.25604      34 2.840029   3.425263   FALSE
## 
## attr(,"class")
## [1] "gofOutlier"
write.csv(acg_dat, "/Users/meaghan/OneDrive - University of Calgary/Preschool_data/MRS/data/acg_data_molal_outliers_removed_n113.csv", row.names = FALSE)

LAG

plot_age <- ggplot(lag_dat, aes(mri_age_y))
plot_age + geom_density()

#age skewed

plot_naa <- ggplot(lag_dat, aes(NAA_conc_molal))
plot_naa + geom_density()

plot_naa + geom_boxplot()

rosnerTest(lag_dat$NAA_conc_molal, k = 10) #outliers val. 19.207, obs. 296 & Val. 10.148, obs. 306 - some issues in spectra, exclude from analysis for NAA
## $distribution
## [1] "Normal"
## 
## $statistic
##      R.1      R.2      R.3      R.4      R.5      R.6      R.7      R.8 
## 4.467031 4.044357 3.599759 3.613937 3.039885 2.720506 2.647000 2.651977 
##      R.9     R.10 
## 2.660497 2.662447 
## 
## $sample.size
## [1] 321
## 
## $parameters
##  k 
## 10 
## 
## $alpha
## [1] 0.05
## 
## $crit.value
##  lambda.1  lambda.2  lambda.3  lambda.4  lambda.5  lambda.6  lambda.7  lambda.8 
##  3.742579  3.741706  3.740829  3.739949  3.739067  3.738181  3.737291  3.736399 
##  lambda.9 lambda.10 
##  3.735503  3.734604 
## 
## $n.outliers
## [1] 2
## 
## $alternative
## [1] "Up to 10 observations are not\n                                 from the same Distribution."
## 
## $method
## [1] "Rosner's Test for Outliers"
## 
## $data
##   [1] 14.16473 13.06076 14.35969 13.27563 16.92063 13.92045 14.90883 14.96705
##   [9] 13.87538 13.79928 15.43119 13.74363 13.42774 14.32941 15.29816 15.36852
##  [17] 14.96343 14.01462 15.11779 16.29892 15.28430 15.17137 15.34158 15.28575
##  [25] 13.61378 16.76759 13.63702 14.41263 13.69061 15.64950 14.60942 15.25717
##  [33] 15.68290 14.24180 15.96102 13.19084 14.22355 14.45688 13.89590 16.68967
##  [41] 15.14115 13.88211 14.56337 14.77687 14.00542 14.12098 13.36743 13.70315
##  [49] 14.59050 15.55276 13.33791 14.13365 12.08635 14.08027 13.90090 13.63249
##  [57] 15.72428 11.86871 14.17563 11.96658 12.37485 14.57472 13.54554 14.35760
##  [65] 13.86325 16.89273 14.68644 16.67785 15.93255 16.85334 14.93659 14.21829
##  [73] 13.99673 14.23697 13.86957 13.86406 14.36716 14.41372 12.45632 14.76312
##  [81] 14.75832 14.92673 13.88514 14.24089 15.42366 14.86462 14.32717 14.96667
##  [89] 13.72426 14.16948 15.16894 13.93828 12.92563 13.94977 13.64160 14.55102
##  [97] 14.51660 12.52392 14.21054 13.93670 16.19979 14.45536 14.10612 13.35969
## [105] 16.58862 14.41649 14.80811 14.26277 12.37256 15.44893 13.43728 13.57859
## [113] 13.72929 14.86513 11.95501 14.39891 14.11323 14.14168 14.73091 14.40429
## [121] 13.19741 13.71443 14.26528 14.92876 14.67063 14.35362 15.12298 13.95432
## [129] 13.68703 13.44733 14.01869 14.75568 14.74986 11.38458 12.08162 11.74260
## [137] 13.93366 14.16181 14.38636 14.40623 14.15278 14.17590 15.46643 14.09616
## [145] 13.62893 14.14878 13.54441 13.09883 12.79401 14.64345 14.02714 14.59131
## [153] 14.51512 15.41640 13.41648 13.54853 13.78080 13.79989 13.51420 14.77851
## [161] 15.60295 13.13390 13.99763 14.27699 15.00866 15.01634 12.09852 14.98452
## [169] 13.60870 16.34165 11.91384 12.33377 14.09992 15.06924 13.87378 13.47019
## [177] 14.02295 18.06089 14.25601 14.37746 13.24330 15.14502 13.95585 13.43714
## [185] 14.30145 16.30323 13.44700 14.70962 13.78765 13.43206 14.86922 14.73623
## [193] 14.54106 13.30455 14.61897 13.91250 15.44818 13.00756 14.21959 16.58017
## [201] 14.30707 15.14338 12.83493 13.96649 15.41109 13.77502 12.24289 15.19071
## [209] 13.32940 13.11120 15.00722 15.15742 13.96143 14.70099 14.27321 15.14983
## [217] 14.79365 15.19557 14.84745 13.19832 15.12751 13.95088 13.28456 14.90679
## [225] 12.94820 13.28282 14.39416 13.76430 14.91352 15.55618 15.20787 14.89958
## [233] 14.31622 14.53646 15.29391 13.51980 15.05685 14.73951 15.01094 15.61748
## [241] 14.09588 14.63147 15.59956 15.33848 14.07436 14.30913 13.94217 16.47833
## [249] 15.29907 13.85813 15.73137 15.08295 16.15169 14.28221 15.30824 14.47596
## [257] 14.47137 14.28477 15.33016 14.75400 14.54703 14.66072 13.49926 14.47973
## [265] 13.83620 14.83537 15.05792 14.18100 13.68113 13.28288 14.53224 13.92267
## [273] 14.80805 13.25470 14.60733 13.77756 15.36693 16.13917 14.36568 14.54767
## [281] 14.51644 16.31083 15.95003 13.41067 15.50889 13.79439 14.53462 14.83200
## [289] 14.42306 13.93431 13.79584 14.91428 12.66542 14.98849 15.38162 19.20747
## [297] 13.74790 17.99338 14.95991 13.93438 15.42154 13.27326 13.14415 14.47594
## [305] 13.99067 10.14812 13.73230 14.91913 15.85376 13.23125 13.40927 14.90803
## [313] 14.30314 14.54287 14.33013 16.55184 14.94880 14.30460 15.31501 14.46473
## [321] 16.31675
## 
## $data.name
## [1] "lag_dat$NAA_conc_molal"
## 
## $bad.obs
## [1] 0
## 
## $all.stats
##    i   Mean.i      SD.i    Value Obs.Num    R.i+1 lambda.i+1 Outlier
## 1  0 14.39149 1.0781156 19.20747     296 4.467031   3.742579    TRUE
## 2  1 14.37644 1.0454867 10.14812     306 4.044357   3.741706    TRUE
## 3  2 14.38970 1.0198436 18.06089     178 3.599759   3.740829   FALSE
## 4  3 14.37815 1.0003560 17.99338     298 3.613937   3.739949   FALSE
## 5  4 14.36675 0.9810139 11.38458     134 3.039885   3.739067   FALSE
## 6  5 14.37619 0.9680502 11.74260     136 2.720506   3.738181   FALSE
## 7  6 14.38455 0.9580956 16.92063       5 2.647000   3.737291   FALSE
## 8  7 14.37647 0.9488234 16.89273      66 2.651977   3.736399   FALSE
## 9  8 14.36843 0.9395707 11.86871      58 2.660497   3.735503   FALSE
## 10 9 14.37644 0.9303092 16.85334      70 2.662447   3.734604   FALSE
## 
## attr(,"class")
## [1] "gofOutlier"
lag_dat["296", "NAA_conc_molal"]<- NA
lag_dat["306", "NAA_conc_molal"]<- NA

plot_cre <- ggplot(lag_dat, aes(Cre_PCr_conc_molal))
plot_cre + geom_density()

plot_cre + geom_boxplot()

rosnerTest(lag_dat$Cre_PCr_conc_molal, k = 10) #outliers val. 14.829, obs. 241 (PS17_047) - structure in resid. & val. 14.09, obs. 178 (PS17_045, spectrum looked fine), excluding only PS17_047 for now
## $distribution
## [1] "Normal"
## 
## $statistic
##      R.1      R.2      R.3      R.4      R.5      R.6      R.7      R.8 
## 4.835902 4.220654 3.332561 3.278883 3.286432 3.307998 2.911341 2.895521 
##      R.9     R.10 
## 2.767415 2.775535 
## 
## $sample.size
## [1] 321
## 
## $parameters
##  k 
## 10 
## 
## $alpha
## [1] 0.05
## 
## $crit.value
##  lambda.1  lambda.2  lambda.3  lambda.4  lambda.5  lambda.6  lambda.7  lambda.8 
##  3.742579  3.741706  3.740829  3.739949  3.739067  3.738181  3.737291  3.736399 
##  lambda.9 lambda.10 
##  3.735503  3.734604 
## 
## $n.outliers
## [1] 2
## 
## $alternative
## [1] "Up to 10 observations are not\n                                 from the same Distribution."
## 
## $method
## [1] "Rosner's Test for Outliers"
## 
## $data
##   [1]  9.979902  8.752796  9.874363  8.961297 11.510180 10.121568 10.418593
##   [8] 11.028922 10.171413  9.929791 10.871830 11.064228 10.111309  9.319128
##  [15] 10.662947 10.890644 10.197139  9.919287 10.426231 11.325982 10.557574
##  [22] 11.088082 11.039646 11.172931  8.936929 12.109674 10.000368 10.420059
##  [29]  9.572731 10.511121  9.984628 10.695378 11.039338  9.822119 11.097416
##  [36]  9.239262  9.742404 10.192047 10.845538 12.117826 10.389433  9.535353
##  [43] 10.010295 10.823639  9.876324 10.019161  8.608303  9.442680 11.444015
##  [50] 11.419644  9.408142 10.139106  8.483950  9.792171 10.032519  8.968350
##  [57] 11.562417  9.156420 10.285179  7.449145 10.159106 10.623228 10.245682
##  [64] 11.514351 10.182653 11.786704 11.560836 12.034078 11.166365 11.463697
##  [71] 10.397742 10.176308 11.796639 10.733740 10.174742 10.465407  9.824296
##  [78] 11.420357  7.832529 10.284534 12.513826 11.576531 10.541650 10.808215
##  [85] 10.483803 10.041962 10.382300  9.976776 10.765842 10.520082 11.422196
##  [92] 10.843472  9.968723 10.781936  9.841352  9.930816 11.494602  9.908886
##  [99] 10.765624 10.122441 11.914466 10.403873 10.056177  9.179081 11.252406
## [106] 10.011672 10.666794  9.981113  9.225683 10.367639 10.879950 11.335061
## [113]  9.425501 10.725975 10.241112  9.317008  9.889557  9.560765 10.593153
## [120] 11.263291  8.935236 10.137273 10.048855 11.211108 10.481657 10.028783
## [127] 11.460008  9.844337 10.207849  8.876822  9.557376  9.418484  9.631575
## [134] 11.004550  7.414051  8.984416  9.798072 10.236338 10.150536 10.108983
## [141] 10.653763 10.137920 10.541536  9.525275  9.886573  9.896683  9.185657
## [148]  9.642602  7.281217  9.968306 10.715037 10.739283 11.039381 11.250002
## [155]  9.759444  9.643949  9.675426  9.621214  9.141950 10.356422 10.788507
## [162]  9.362209  9.822128  9.645301 10.157630 10.381588  8.913720 11.109833
## [169]  9.890247 12.488886  9.115647  9.814255 10.627685 11.734993 10.295594
## [176] 10.173568  9.988731 14.090723 10.829617 10.247883 10.533434 10.999808
## [183] 10.001382  9.229915 10.354994 11.296388  9.994843 10.904378  9.345731
## [190]  9.618365  9.759703 10.530182 10.844392  8.895843 10.110129 11.423088
## [197] 11.096328  9.958072  9.931007 10.921273  9.154313 10.482202 10.423407
## [204] 10.629254 12.212044  9.913344  8.305599 10.446897  9.841379  9.379922
## [211]  9.844787 10.968766 10.801476  9.891839 10.632212  9.671280 10.872873
## [218] 10.666877 10.645958  9.493899 13.103099  9.512697 10.276565  9.846960
## [225]  9.773848  9.929345 11.034653 10.095542 10.422810 10.640722 10.407355
## [232] 11.450386 10.284795 10.179735 11.263605 10.183927 10.257883 10.937157
## [239] 10.529431 11.012918 14.829071  9.133491 11.531961 10.173887  9.610385
## [246]  9.349643  9.324653 10.926724 11.044735  8.859819 11.228731 10.816240
## [253] 11.271292  9.861827 10.852278 10.764773  9.380198  9.481786 10.781017
## [260]  8.985104  9.888393  9.434081  8.870117  9.122209  9.542512  9.889146
## [267] 10.769356  9.616033  8.966342  9.255494 10.764599  9.897261  9.427714
## [274]  8.163350 10.017130  9.265501 10.944724 12.004282  9.764853  9.856451
## [281] 10.020483 12.009948 10.044769  9.667192 11.197137  9.964094 10.035066
## [288] 10.812708  9.728669  9.918824  8.843802 10.202045  9.183452  9.822267
## [295] 10.790509 12.454846  8.641534  9.507122 10.361671  9.646740 11.011202
## [302]  9.075410 10.166484  9.737467  9.615290  7.882414 10.377305 10.224366
## [309] 12.449600 10.083783  9.260131 11.009901  9.914846  9.659798 10.610323
## [316] 11.561746 11.077058  9.291690 10.160005  9.934510 11.934339
## 
## $data.name
## [1] "lag_dat$Cre_PCr_conc_molal"
## 
## $bad.obs
## [1] 0
## 
## $all.stats
##    i   Mean.i      SD.i     Value Obs.Num    R.i+1 lambda.i+1 Outlier
## 1  0 10.26105 0.9446048 14.829071     241 4.835902   3.742579    TRUE
## 2  1 10.24678 0.9107457 14.090723     178 4.220654   3.741706    TRUE
## 3  2 10.23473 0.8862592  7.281217     149 3.332561   3.740829   FALSE
## 4  3 10.24402 0.8719682 13.103099     221 3.278883   3.739949   FALSE
## 5  4 10.23500 0.8583617  7.414051     135 3.286432   3.739067   FALSE
## 6  5 10.24393 0.8448557  7.449145      60 3.307998   3.738181   FALSE
## 7  6 10.25280 0.8313244  7.832529      79 2.911341   3.737291   FALSE
## 8  7 10.26051 0.8213002  7.882414     306 2.895521   3.736399   FALSE
## 9  8 10.26810 0.8114875 12.513826      81 2.767415   3.735503   FALSE
## 10 9 10.26091 0.8027211 12.488886     170 2.775535   3.734604   FALSE
## 
## attr(,"class")
## [1] "gofOutlier"
lag_dat["241", "Cre_PCr_conc_molal"]<-NA  #replace outlier value with NA to exclude from analysis

plot_cho <- ggplot(lag_dat, aes(Cho_GPC_PCh_conc_molal))
plot_cho + geom_density()

plot_cho + geom_boxplot()

rosnerTest(lag_dat$Cho_GPC_PCh_conc_molal, k = 8)
## $distribution
## [1] "Normal"
## 
## $statistic
##      R.1      R.2      R.3      R.4      R.5      R.6      R.7      R.8 
## 3.694967 3.121637 3.061660 3.062209 2.979166 3.001536 2.871646 2.733239 
## 
## $sample.size
## [1] 321
## 
## $parameters
## k 
## 8 
## 
## $alpha
## [1] 0.05
## 
## $crit.value
## lambda.1 lambda.2 lambda.3 lambda.4 lambda.5 lambda.6 lambda.7 lambda.8 
## 3.742579 3.741706 3.740829 3.739949 3.739067 3.738181 3.737291 3.736399 
## 
## $n.outliers
## [1] 0
## 
## $alternative
## [1] "Up to 8 observations are not\n                                 from the same Distribution."
## 
## $method
## [1] "Rosner's Test for Outliers"
## 
## $data
##   [1] 2.068181 2.062467 2.386391 2.355035 2.498191 2.256110 2.077815 2.425292
##   [9] 2.398770 1.895982 1.865260 2.184158 2.004899 2.128662 2.236830 2.378977
##  [17] 2.355458 2.266506 2.404117 2.307009 2.052566 2.364941 2.267604 1.940308
##  [25] 2.587651 2.716570 2.469831 2.547177 2.087747 2.170207 2.232363 2.241430
##  [33] 2.441839 2.264953 2.333662 2.359801 2.270399 2.241648 1.920730 2.467142
##  [41] 2.242001 2.227010 2.229581 2.053798 1.939191 2.364334 2.678152 2.362113
##  [49] 1.676867 2.166236 2.306276 2.211050 2.178445 2.361600 2.344510 2.268087
##  [57] 2.667157 2.059386 2.034573 2.045022 2.214431 2.995506 1.846508 2.183946
##  [65] 2.072392 2.429630 2.330074 2.491408 2.265040 2.154392 2.411001 2.046897
##  [73] 2.534412 2.027012 2.376839 2.222266 2.306644 2.292419 2.108096 2.172672
##  [81] 1.839589 1.859195 1.443755 1.640698 2.420001 2.337525 2.238790 2.240345
##  [89] 2.290795 2.399979 2.350376 2.151479 2.142910 2.110996 1.985190 1.928934
##  [97] 2.304258 1.765061 2.137595 1.985591 2.781876 2.123475 1.756957 1.700318
## [105] 1.975316 1.974756 1.999963 2.031005 1.687589 1.708882 2.296279 2.182991
## [113] 2.276755 2.101586 1.837847 2.371711 2.065515 2.499424 2.363551 2.223515
## [121] 2.232809 2.045991 1.947255 2.339091 2.193248 2.528316 2.691920 2.777161
## [129] 2.451074 2.763753 2.377191 2.160428 2.215103 2.214312 2.533792 2.384939
## [137] 2.179386 2.552211 2.673232 2.476304 1.625876 1.926042 2.202001 2.204096
## [145] 2.344069 2.152159 1.895155 1.710908 1.401447 1.946252 2.047116 2.456831
## [153] 2.086612 2.201666 2.190181 2.294616 2.581968 2.225911 2.498999 2.000618
## [161] 2.163682 1.960245 2.204252 2.159811 2.097381 2.254865 1.844443 1.598729
## [169] 1.951488 3.208685 1.836756 1.950559 1.992121 2.281328 2.008449 2.599790
## [177] 2.317712 3.033330 2.698335 2.332814 2.583764 2.234143 2.191218 2.264957
## [185] 1.977277 2.067888 2.356329 2.431505 2.779803 2.150897 2.566469 2.101205
## [193] 2.159059 2.164490 2.152499 2.489624 2.791537 1.673365 2.193421 2.761471
## [201] 1.959300 2.022674 1.970196 2.366537 2.774857 2.129388 1.836375 2.360163
## [209] 2.476086 2.665161 2.258604 2.298965 2.012883 2.201029 2.158646 2.353116
## [217] 1.946114 1.912304 1.692442 2.307008 2.449170 2.132806 2.045527 2.155850
## [225] 2.166616 2.440956 1.783126 2.272129 2.595674 2.284913 2.155251 2.400964
## [233] 2.048543 2.045160 2.356748 2.222524 2.346513 2.378440 2.281038 1.906127
## [241] 1.725583 2.479767 2.262261 2.082256 2.047599 2.138599 2.252565 2.262486
## [249] 1.881557 2.322610 1.489680 2.225178 2.118903 2.390375 2.437536 2.400945
## [257] 1.832693 2.290250 2.007843 2.346441 2.089229 2.428740 2.391081 2.732220
## [265] 1.971298 2.235255 2.116196 1.941129 2.009743 1.999019 1.927678 1.994412
## [273] 2.372282 2.491822 1.886355 2.109659 2.028203 1.807961 2.439615 2.152157
## [281] 1.892138 2.163111 2.376896 1.914835 2.734225 1.944843 1.756721 2.135482
## [289] 2.101386 1.768588 1.974516 1.975407 2.167757 2.112266 2.510332 2.521732
## [297] 2.431779 2.959724 2.046249 1.986009 2.300110 2.368720 2.151772 2.070465
## [305] 2.457353 2.243696 2.105342 2.108232 2.180180 2.349512 1.911659 2.150353
## [313] 2.011103 2.243254 2.417149 1.800195 2.190437 2.875345 2.116343 2.319437
## [321] 2.422916
## 
## $data.name
## [1] "lag_dat$Cho_GPC_PCh_conc_molal"
## 
## $bad.obs
## [1] 0
## 
## $all.stats
##   i   Mean.i      SD.i    Value Obs.Num    R.i+1 lambda.i+1 Outlier
## 1 0 2.208213 0.2707661 3.208685     170 3.694967   3.742579   FALSE
## 2 1 2.205087 0.2653235 3.033330     178 3.121637   3.741706   FALSE
## 3 2 2.202490 0.2616371 1.401447     149 3.061660   3.740829   FALSE
## 4 3 2.205009 0.2581460 2.995506      62 3.062209   3.739949   FALSE
## 5 4 2.202516 0.2546890 1.443755      83 2.979166   3.739067   FALSE
## 6 5 2.204917 0.2514736 2.959724     298 3.001536   3.738181   FALSE
## 7 6 2.202521 0.2482341 1.489680     251 2.871646   3.737291   FALSE
## 8 7 2.204791 0.2453333 2.875345     318 2.733239   3.736399   FALSE
## 
## attr(,"class")
## [1] "gofOutlier"
plot_ins <- ggplot(lag_dat, aes(Ins_conc_molal))
plot_ins + geom_density()

plot_ins + geom_boxplot()

rosnerTest(lag_dat$Ins_conc_molal, k = 8) #outliers: val. 14.077, obs. 241 (PS17_047) - structure in resid; val. 11.777, obs. 221 (PS16_014) - spectrum fine, keep for now; val. 1.542, obs. 197 (PS18_028) - noisy spectrum, spiky residuals, exclude; val. 2.04, obs. 62 (PS14_111) - spectrum fine, keep for now
## $distribution
## [1] "Normal"
## 
## $statistic
##      R.1      R.2      R.3      R.4      R.5      R.6      R.7      R.8 
## 5.981441 4.463363 4.077869 3.761599 3.384729 3.303295 3.148162 3.105675 
## 
## $sample.size
## [1] 321
## 
## $parameters
## k 
## 8 
## 
## $alpha
## [1] 0.05
## 
## $crit.value
## lambda.1 lambda.2 lambda.3 lambda.4 lambda.5 lambda.6 lambda.7 lambda.8 
## 3.742579 3.741706 3.740829 3.739949 3.739067 3.738181 3.737291 3.736399 
## 
## $n.outliers
## [1] 4
## 
## $alternative
## [1] "Up to 8 observations are not\n                                 from the same Distribution."
## 
## $method
## [1] "Rosner's Test for Outliers"
## 
## $data
##   [1]  8.084411  7.069186  6.913031  6.915494  7.726417  8.046174  6.937152
##   [8]  4.793633  5.952276  6.984233  6.400936  4.413732  7.176021  7.448535
##  [15]  7.547442  6.541306  8.400250  5.957678  6.965664  7.688786  6.245160
##  [22]  7.360921  7.925104  4.710055  6.035369  7.937845  7.166536  7.269577
##  [29]  5.736145  7.751891  6.217281  7.121735  7.138640  6.104017  7.736566
##  [36]  6.539681  4.958196  7.275083  6.485777  6.598573  7.153802  7.352086
##  [43]  7.119035  6.146843  6.038556  6.638197  4.629802  6.187281  6.553798
##  [50]  5.945866  7.474194  5.559932  5.743412  6.163849  8.347824  4.432658
##  [57]  8.169688  5.045280  8.873237  6.213139  6.400722  2.040275  8.158702
##  [64]  6.760986  6.468052  7.273556  6.742605  7.552488  7.843304  5.876394
##  [71]  8.162408  4.963944  7.945325  5.440451  4.854219  4.834565  5.839066
##  [78]  5.839637  5.890621  5.871859  6.963936  5.579543  5.190653  7.062975
##  [85]  6.065434  7.046003  6.865074  6.210505  6.024185  6.121049  8.133551
##  [92]  6.863480  6.537634  7.457231  5.008241  6.892333  5.936366  3.826364
##  [99]  7.075620  5.567096  5.248858  5.491554  4.529573  7.011630  5.896171
## [106]  6.268198  5.292489  5.663161  6.753982  5.643647  6.788746  5.437125
## [113]  3.078293  7.511207  2.568354  8.286250  5.309958  6.053735  8.092627
## [120]  7.122038  6.509101  6.564424  5.780997  6.770400  5.727577  7.829891
## [127]  6.206182  7.258719  6.558572  4.787960  6.132452  6.029208  6.663446
## [134]  7.630692  4.608510  5.731585  5.435869  5.209567  2.734199  5.857696
## [141]  5.638088  6.135153  6.888234  5.236533  5.655180  4.948126  4.992038
## [148]  6.025793  2.972647  6.870427  7.728083  7.800580  7.334580  6.676194
## [155]  6.170230  5.196197  7.073851  6.127912  7.232895  4.960469  6.612268
## [162]  6.088183  5.798318  6.543746  6.997835  6.987522  4.250557  4.637685
## [169]  6.778603  7.502499  4.679992  5.270903  6.367636  6.301188  7.918660
## [176]  6.533953  7.497948  8.801777  6.797100  5.264895  5.014240  6.141243
## [183]  5.472028  3.885218  4.816043  6.248540  4.626087  6.007857  6.538200
## [190]  7.586651  7.411192  7.320786  7.682231  7.120925  6.686811  7.373729
## [197]  1.541994  6.314608  6.840831  7.095049  6.617076  6.056406  5.882616
## [204]  4.775951  5.731427  5.979244  6.418407  6.735628  5.243289  5.869558
## [211]  5.226573  6.275003  7.093997  7.541310  5.731832  4.649731  6.150065
## [218]  6.428658  7.426624  6.729085 11.777066  5.591387  3.080430  8.027587
## [225]  7.066068  6.661947  6.600256  6.955873  6.949031  6.824263  7.605754
## [232]  6.414457  6.094616  6.805794  5.865261  8.270934  4.935241  6.818562
## [239]  5.441384  7.674150 14.076805  4.939637  6.137708  5.300423  6.306576
## [246]  6.245176  7.663984  7.103391  8.335725  3.649002  9.145931  7.122485
## [253]  8.331023  7.199990  8.067478  7.557525  5.521439  6.036198  6.082461
## [260]  7.524328  5.926927  5.705055  5.399996  6.232379  5.707127  4.605042
## [267]  4.333050  5.607343  5.940738  5.740379  6.177392  6.685727  5.210419
## [274]  6.865329  6.749236  5.464544  6.496419  6.705573  7.274507  6.221882
## [281]  5.994352  6.059500  7.660503  6.660801  7.439851  6.777698  5.282780
## [288]  8.773125  6.173129  5.673978  5.639679  5.174471  8.615918  5.774767
## [295]  5.411721  7.832813  7.925781  4.996084  6.276000  6.161567  6.484734
## [302]  5.594141  6.793109  6.401047  9.077790  8.562826  5.674753  6.368745
## [309]  8.425749  6.658481  4.433254  5.482190  7.149805  6.079902  6.139307
## [316]  7.111624  5.971685  6.120964  6.571466  5.418973  6.379001
## 
## $data.name
## [1] "lag_dat$Ins_conc_molal"
## 
## $bad.obs
## [1] 0
## 
## $all.stats
##   i   Mean.i     SD.i     Value Obs.Num    R.i+1 lambda.i+1 Outlier
## 1 0 6.383077 1.286267 14.076805     241 5.981441   3.742579    TRUE
## 2 1 6.359034 1.213890 11.777066     221 4.463363   3.741706    TRUE
## 3 2 6.342050 1.177099  1.541994     197 4.077869   3.740829    TRUE
## 4 3 6.357144 1.147615  2.040275      62 3.761599   3.739949    TRUE
## 5 4 6.370762 1.123401  2.568354     115 3.384729   3.739067   FALSE
## 6 5 6.382795 1.104533  2.734199     139 3.303295   3.738181   FALSE
## 7 6 6.394378 1.086898  2.972647     149 3.148162   3.737291   FALSE
## 8 7 6.405275 1.071259  3.078293     113 3.105675   3.736399   FALSE
## 
## attr(,"class")
## [1] "gofOutlier"
lag_dat["241", "Ins_conc_molal"] <- NA #replace outlier value with NA to exclude from analysis
lag_dat["197", "Ins_conc_molal"] <- NA

plot_glx <- ggplot(lag_dat, aes(Glu_Gln_conc_molal))
plot_glx + geom_density()

plot_glx + geom_boxplot()

rosnerTest(lag_dat$Glu_Gln_conc_molal, k = 7)#outliers: val. 10.09, obs. 241 (PS17_047) - structure in resid; val. 31.2899, obs. 296 (PS21_030) - big spike in resid. around 2.1, exclude.
## $distribution
## [1] "Normal"
## 
## $statistic
##      R.1      R.2      R.3      R.4      R.5      R.6      R.7 
## 4.954047 4.247538 3.396313 3.439741 3.463175 3.267721 3.232088 
## 
## $sample.size
## [1] 321
## 
## $parameters
## k 
## 7 
## 
## $alpha
## [1] 0.05
## 
## $crit.value
## lambda.1 lambda.2 lambda.3 lambda.4 lambda.5 lambda.6 lambda.7 
## 3.742579 3.741706 3.740829 3.739949 3.739067 3.738181 3.737291 
## 
## $n.outliers
## [1] 2
## 
## $alternative
## [1] "Up to 7 observations are not\n                                 from the same Distribution."
## 
## $method
## [1] "Rosner's Test for Outliers"
## 
## $data
##   [1] 23.58311 20.72833 21.43465 20.79838 24.55542 19.38008 22.13533 21.98148
##   [9] 18.58332 21.80094 23.39061 22.39400 22.25994 19.84839 22.62240 19.74196
##  [17] 19.55718 20.18644 20.93614 24.01329 22.14911 23.96545 23.49753 23.68093
##  [25] 20.83139 25.55075 18.83744 23.02848 20.10324 24.84118 23.11824 25.62522
##  [33] 25.14532 20.60338 23.10183 18.68727 19.90887 22.63726 23.44172 23.97287
##  [41] 20.79684 21.03860 21.41871 22.38334 19.07457 19.80466 21.54623 21.54133
##  [49] 22.80842 25.04145 20.65555 22.02043 18.52301 21.98565 22.88652 20.99602
##  [57] 23.64194 20.45706 24.44658 22.40116 19.22839 24.77315 20.40702 22.32472
##  [65] 22.79321 23.27919 23.66026 23.59194 23.28915 23.15967 20.72270 22.14974
##  [73] 23.87191 27.09158 23.12072 25.11318 22.08527 21.17324 14.25503 21.44083
##  [81] 23.08444 21.47161 23.79592 26.39922 24.52352 23.18206 22.57951 22.53444
##  [89] 24.24539 23.38651 21.39081 23.50524 21.84539 23.16617 20.66627 20.68771
##  [97] 21.33913 23.82111 24.17549 20.43359 24.17777 21.80255 21.79162 21.34246
## [105] 26.43216 21.02534 21.61161 22.57479 21.58113 22.47585 21.89183 22.86451
## [113] 21.63056 23.83333 23.16833 22.80824 21.74170 19.79860 22.68433 25.43932
## [121] 20.46600 19.55287 21.71921 23.91575 19.94617 22.54020 25.12373 20.26601
## [129] 20.42848 18.40115 22.31740 21.34145 21.18639 24.46723 17.55526 19.17975
## [137] 22.14220 21.22080 20.62138 19.29726 22.33080 22.36491 23.64798 24.04583
## [145] 25.02751 22.02316 20.84404 22.82850 23.75262 19.93399 21.93649 21.46368
## [153] 23.51970 23.79151 22.29261 22.55887 19.88319 23.52534 18.26039 22.44651
## [161] 23.87196 21.33171 21.19221 21.21896 25.13735 21.72890 21.00404 24.38403
## [169] 22.77448 24.07955 20.92792 22.79597 21.53526 25.58587 20.23221 20.38536
## [177] 19.64011 29.12643 21.85437 22.81675 24.27197 20.09271 21.28381 21.20096
## [185] 20.61092 26.55674 20.24095 22.22341 21.17598 21.05662 20.70603 22.17497
## [193] 22.22738 17.77073 21.62541 23.02400 23.16615 17.38515 20.64300 21.97048
## [201] 23.41849 22.42124 20.67899 21.35918 23.92414 18.75710 19.72431 22.35249
## [209] 23.13971 19.66645 21.97753 21.86935 22.15373 18.04793 20.63749 19.81250
## [217] 22.91987 21.99328 21.86022 19.23943 22.18618 21.97414 22.10059 21.96999
## [225] 21.38803 21.08812 21.12026 20.43893 20.01495 22.48505 23.44854 23.27656
## [233] 22.72471 21.03119 22.44160 19.23947 24.24070 21.72010 21.19354 21.24039
## [241] 10.09497 19.32810 23.68178 22.11865 18.91616 19.59606 19.02554 22.92979
## [249] 26.99382 18.91880 22.36158 21.71252 22.51005 20.90044 23.69281 22.60270
## [257] 20.69658 19.00693 19.97580 15.08181 17.70438 17.08767 17.40514 16.78328
## [265] 20.40818 19.57305 21.57031 21.73287 18.10301 16.64465 21.45075 20.54401
## [273] 19.32920 14.36555 21.66594 17.20659 21.57909 24.57946 17.49328 18.35123
## [281] 21.22156 24.11680 22.51873 19.19853 23.08553 19.03798 23.38476 21.59376
## [289] 19.09657 21.28636 18.16122 19.70762 18.41625 18.73419 22.65583 31.28994
## [297] 18.89292 17.32671 18.83142 19.60466 24.66522 18.35557 22.92964 18.49429
## [305] 21.56271 18.98307 20.17174 24.36413 25.12361 21.76663 19.89145 21.01574
## [313] 20.06260 25.41609 26.15708 23.69646 23.29001 20.54746 23.65148 20.44966
## [321] 28.51150
## 
## $data.name
## [1] "lag_dat$Glu_Gln_conc_molal"
## 
## $bad.obs
## [1] 0
## 
## $all.stats
##   i   Mean.i     SD.i    Value Obs.Num    R.i+1 lambda.i+1 Outlier
## 1 0 21.68857 2.340229 10.09497     241 4.954047   3.742579    TRUE
## 2 1 21.72480 2.251926 31.28994     296 4.247538   3.741706    TRUE
## 3 2 21.69482 2.190549 14.25503      79 3.396313   3.740829   FALSE
## 4 3 21.71821 2.153715 29.12643     178 3.439741   3.739949   FALSE
## 5 4 21.69484 2.116351 14.36555     274 3.463175   3.739067   FALSE
## 6 5 21.71804 2.078963 28.51150     321 3.267721   3.738181   FALSE
## 7 6 21.69647 2.046559 15.08181     260 3.232088   3.737291   FALSE
## 
## attr(,"class")
## [1] "gofOutlier"
lag_dat["241", "Glu_Gln_conc_molal"] <- NA
lag_dat["296", "Glu_Gln_conc_molal"] <- NA

#exclude PS17_047 from all analyses due to outlier status on multiple metabolites and issues in spectrum (noisy, structure in resid)
lag_dat["241", "NAA_conc_molal"] <- NA
lag_dat["241", "Cre_PCr_conc_molal"] <- NA
lag_dat["241", "Cho_GPC_PCh_conc_molal"] <- NA

write.csv(lag_dat, "/Users/meaghan/OneDrive - University of Calgary/Preschool_data/MRS/data/lag_data_molal_outliers_removed_n320.csv", row.names = FALSE)

##Descriptive stats

#ACG - count individual subject N and N of females
acg_sub_sex <- acg_dat %>%
  dplyr::select(subj_id, female)
acg_unique_subs<-unique(acg_sub_sex)
length(acg_unique_subs$subj_id)
## [1] 102
sum(acg_unique_subs$female)
## [1] 47
describe(acg_dat)
##                        vars   n    mean      sd  median trimmed     mad     min
## study_code*               1 113   57.00   32.76   57.00   57.00   41.51    1.00
## subj_id*                  2 113   51.42   28.89   54.00   51.46   35.58    1.00
## spec_id                   3 113 5072.42 3034.48 3244.00 4708.09 1254.28 2058.00
## date_scan*                4 113   49.21   28.71   50.00   49.32   37.06    1.00
## hand*                     5 113    4.54    1.04    5.00    4.78    0.00    1.00
## female                    6 113    0.47    0.50    0.00    0.46    0.00    0.00
## mri_age_y                 7 113    4.06    1.10    3.68    3.97    1.00    2.34
## folder*                   8 113   57.00   32.76   57.00   57.00   41.51    1.00
## File_ID*                  9 113   57.00   32.76   57.00   57.00   41.51    1.00
## fGM                      10 113    0.77    0.05    0.78    0.78    0.05    0.55
## fWM                      11 113    0.04    0.03    0.03    0.03    0.02    0.00
## fCSF                     12 113    0.19    0.05    0.18    0.19    0.04    0.09
## NAA_conc_molal           13 113   14.89    1.32   15.03   14.94    1.25   11.24
## Cre_PCr_conc_molal       14 113   12.04    1.07   12.04   12.06    0.84    8.41
## Cho_GPC_PCh_conc_molal   15 112    2.50    0.34    2.48    2.49    0.30    1.57
## Ins_conc_molal           16 112    7.64    1.53    7.71    7.72    1.46    2.99
## Glu_conc_molal           17 113   22.38    2.25   22.50   22.46    2.01   16.03
## Glu_Gln_conc_molal       18 113   26.83    2.62   26.72   26.81    2.64   20.84
##                             max   range  skew kurtosis     se
## study_code*              113.00  112.00  0.00    -1.23   3.08
## subj_id*                 102.00  101.00 -0.05    -1.15   2.72
## spec_id                11448.00 9390.00  0.87    -0.81 285.46
## date_scan*                97.00   96.00 -0.04    -1.29   2.70
## hand*                      6.00    5.00 -2.07     3.38   0.10
## female                     1.00    1.00  0.12    -2.00   0.05
## mri_age_y                  7.36    5.03  0.72    -0.49   0.10
## folder*                  113.00  112.00  0.00    -1.23   3.08
## File_ID*                 113.00  112.00  0.00    -1.23   3.08
## fGM                        0.87    0.32 -1.14     1.59   0.01
## fWM                        0.29    0.29  4.07    25.28   0.00
## fCSF                       0.36    0.27  0.90     0.92   0.00
## NAA_conc_molal            18.30    7.05 -0.37     0.21   0.12
## Cre_PCr_conc_molal        14.76    6.35 -0.36     0.83   0.10
## Cho_GPC_PCh_conc_molal     3.33    1.76  0.16     0.08   0.03
## Ins_conc_molal            10.89    7.90 -0.51     0.35   0.14
## Glu_conc_molal            28.20   12.17 -0.31     0.40   0.21
## Glu_Gln_conc_molal        34.26   13.42  0.13    -0.01   0.25
#LAG - count individual subject N and N of females
lag_sub_sex <- lag_dat %>%
  dplyr::select(subj_id, female)
lag_unique_subs<-unique(lag_sub_sex)
length(lag_unique_subs$subj_id)
## [1] 95
sum(lag_unique_subs$female)
## [1] 47
describe(lag_dat)
##                        vars   n    mean      sd  median trimmed     mad     min
## study_code*               1 321  160.94   92.72  161.00  161.00  118.61    1.00
## subj_id*                  2 321   40.27   24.55   41.00   39.16   28.17    1.00
## spec_id                   3 320 7744.26 3536.93 6379.50 7565.04 3846.61 3118.00
## date_scan*                4 321  143.73   82.33  143.00  143.92  103.78    1.00
## hand*                     5 321    4.64    0.90    5.00    4.88    0.00    1.00
## female                    6 321    0.47    0.50    0.00    0.46    0.00    0.00
## mri_age_y                 7 321    5.95    1.90    5.39    5.79    1.79    2.41
## folder*                   8 321  161.00   92.81  161.00  161.00  118.61    1.00
## File_ID*                  9 321  161.00   92.81  161.00  161.00  118.61    1.00
## fGM                      10 321    0.62    0.10    0.64    0.63    0.10    0.30
## fWM                      11 321    0.35    0.11    0.34    0.34    0.11    0.00
## fCSF                     12 321    0.03    0.02    0.03    0.03    0.01    0.00
## NAA_conc_molal           13 318   14.39    1.02   14.36   14.37    0.88   11.38
## Cre_PCr_conc_molal       14 320   10.25    0.91   10.17   10.23    0.83    7.28
## Cho_GPC_PCh_conc_molal   15 320    2.21    0.27    2.20    2.20    0.24    1.40
## Ins_conc_molal           16 319    6.37    1.19    6.37    6.39    1.08    2.04
## Glu_conc_molal           17 321   17.88    1.90   17.87   17.88    1.73   10.63
## Glu_Gln_conc_molal       18 319   21.69    2.19   21.74   21.72    2.02   14.26
##                             max    range  skew kurtosis     se
## study_code*              320.00   319.00  0.00    -1.21   5.17
## subj_id*                  95.00    94.00  0.28    -0.83   1.37
## spec_id                14205.00 11087.00  0.35    -1.37 197.72
## date_scan*               284.00   283.00  0.00    -1.19   4.60
## hand*                      6.00     5.00 -2.35     4.48   0.05
## female                     1.00     1.00  0.13    -1.99   0.03
## mri_age_y                 11.13     8.72  0.66    -0.46   0.11
## folder*                  321.00   320.00  0.00    -1.21   5.18
## File_ID*                 321.00   320.00  0.00    -1.21   5.18
## fGM                        0.80     0.50 -0.75     0.40   0.01
## fWM                        0.66     0.65  0.43     0.20   0.01
## fCSF                       0.21     0.21  3.12    16.87   0.00
## NAA_conc_molal            18.06     6.68  0.20     0.95   0.06
## Cre_PCr_conc_molal        14.09     6.81  0.17     1.33   0.05
## Cho_GPC_PCh_conc_molal     3.21     1.81  0.22     0.80   0.02
## Ins_conc_molal            11.78     9.74 -0.13     1.69   0.07
## Glu_conc_molal            24.87    14.24 -0.07     1.22   0.11
## Glu_Gln_conc_molal        29.13    14.87 -0.12     0.74   0.12
#Calculate total Ns from both datasets combined
all_dat<-bind_rows(acg_dat, lag_dat)
length(unique(all_dat$File_ID))
## [1] 434
length(unique(all_dat$subj_id))
## [1] 127
all_sub_sex <- all_dat %>%
  dplyr::select(subj_id, female)
all_unique_subs<-unique(all_sub_sex)
sum(all_unique_subs$female)
## [1] 62
describe(all_dat)
##                        vars   n    mean      sd  median trimmed     mad     min
## study_code*               1 434  216.94  124.93  216.50  216.97  160.12    1.00
## subj_id*                  2 434   53.87   33.15   54.00   52.20   38.55    1.00
## spec_id                   3 433 7046.99 3606.23 5832.00 6818.88 3964.47 2058.00
## date_scan*                4 434  184.27  105.60  183.00  184.45  133.43    1.00
## hand*                     5 434    4.62    0.93    5.00    4.86    0.00    1.00
## female                    6 434    0.47    0.50    0.00    0.46    0.00    0.00
## mri_age_y                 7 434    5.46    1.92    4.95    5.28    1.68    2.34
## folder*                   8 434  217.04  125.00  217.50  217.05  160.12    1.00
## File_ID*                  9 434  217.50  125.43  217.50  217.50  160.86    1.00
## fGM                      10 434    0.66    0.11    0.67    0.67    0.11    0.30
## fWM                      11 434    0.27    0.17    0.29    0.27    0.16    0.00
## fCSF                     12 434    0.07    0.08    0.03    0.06    0.03    0.00
## NAA_conc_molal           13 431   14.52    1.13   14.48   14.51    0.99   11.24
## Cre_PCr_conc_molal       14 433   10.71    1.24   10.54   10.65    1.10    7.28
## Cho_GPC_PCh_conc_molal   15 432    2.28    0.32    2.26    2.27    0.26    1.40
## Ins_conc_molal           16 431    6.70    1.40    6.66    6.68    1.28    2.04
## Glu_conc_molal           17 434   19.05    2.81   18.53   18.89    2.53   10.63
## Glu_Gln_conc_molal       18 432   23.04    3.23   22.54   22.85    2.77   14.26
##                             max    range  skew kurtosis     se
## study_code*              432.00   431.00  0.00    -1.21   6.00
## subj_id*                 127.00   126.00  0.33    -0.78   1.59
## spec_id                14205.00 12147.00  0.45    -1.19 173.30
## date_scan*               366.00   365.00  0.00    -1.18   5.07
## hand*                      6.00     5.00 -2.29     4.27   0.04
## female                     1.00     1.00  0.13    -1.99   0.02
## mri_age_y                 11.13     8.79  0.79    -0.09   0.09
## folder*                  433.00   432.00  0.00    -1.20   6.00
## File_ID*                 434.00   433.00  0.00    -1.21   6.02
## fGM                        0.87     0.57 -0.60     0.10   0.01
## fWM                        0.66     0.66 -0.09    -0.87   0.01
## fCSF                       0.36     0.36  1.37     0.83   0.00
## NAA_conc_molal            18.30     7.05  0.10     0.57   0.05
## Cre_PCr_conc_molal        14.76     7.48  0.44     0.20   0.06
## Cho_GPC_PCh_conc_molal     3.33     1.93  0.45     0.70   0.02
## Ins_conc_molal            11.78     9.74  0.08     0.70   0.07
## Glu_conc_molal            28.20    17.57  0.49     0.13   0.13
## Glu_Gln_conc_molal        34.26    20.00  0.54     0.34   0.16

LME modelling ACG

run mixed-effects models for ACG with participant as random effect, age, sex and age x sex interaction as fixed effects NAA run on acg_data_molal_outliers_removed_n113.csv

acg_naa_age<- lmer(NAA_conc_molal ~ mri_age_y + (1 | subj_id), data=acg_dat, REML=FALSE)
summary(acg_naa_age)
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
##   method [lmerModLmerTest]
## Formula: NAA_conc_molal ~ mri_age_y + (1 | subj_id)
##    Data: acg_dat
## 
##      AIC      BIC   logLik deviance df.resid 
##    383.3    394.2   -187.6    375.3      109 
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -1.53495 -0.39238  0.05273  0.38044  1.68117 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  subj_id  (Intercept) 1.1187   1.0577  
##  Residual             0.6121   0.7824  
## Number of obs: 113, groups:  subj_id, 102
## 
## Fixed effects:
##             Estimate Std. Error       df t value Pr(>|t|)    
## (Intercept)  13.8379     0.4588 108.1511   30.16   <2e-16 ***
## mri_age_y     0.2513     0.1079 105.1446    2.33   0.0217 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##           (Intr)
## mri_age_y -0.959
acg_naa_age_sex<- lmer(NAA_conc_molal ~ mri_age_y + female + (1 | subj_id), data=acg_dat, REML=FALSE)
summary(acg_naa_age_sex)
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
##   method [lmerModLmerTest]
## Formula: NAA_conc_molal ~ mri_age_y + female + (1 | subj_id)
##    Data: acg_dat
## 
##      AIC      BIC   logLik deviance df.resid 
##    385.2    398.8   -187.6    375.2      108 
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -1.54531 -0.40334  0.05382  0.37656  1.67013 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  subj_id  (Intercept) 1.1189   1.0578  
##  Residual             0.6113   0.7819  
## Number of obs: 113, groups:  subj_id, 102
## 
## Fixed effects:
##              Estimate Std. Error        df t value Pr(>|t|)    
## (Intercept)  13.85437    0.46506 109.75797  29.791   <2e-16 ***
## mri_age_y     0.25366    0.10837 104.64089   2.341   0.0211 *  
## female       -0.05648    0.26075  99.22276  -0.217   0.8289    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##           (Intr) mr_g_y
## mri_age_y -0.926       
## female    -0.165 -0.098
anova(acg_naa_age,acg_naa_age_sex)
## Data: acg_dat
## Models:
## acg_naa_age: NAA_conc_molal ~ mri_age_y + (1 | subj_id)
## acg_naa_age_sex: NAA_conc_molal ~ mri_age_y + female + (1 | subj_id)
##                 npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)
## acg_naa_age        4 383.25 394.16 -187.62   375.25                     
## acg_naa_age_sex    5 385.20 398.84 -187.60   375.20 0.0469  1     0.8285
acg_naa_ageXsex<- lmer(NAA_conc_molal ~ mri_age_y + female + mri_age_y:female + (1 | subj_id), data=acg_dat, REML=FALSE)
summary(acg_naa_ageXsex)
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
##   method [lmerModLmerTest]
## Formula: NAA_conc_molal ~ mri_age_y + female + mri_age_y:female + (1 |  
##     subj_id)
##    Data: acg_dat
## 
##      AIC      BIC   logLik deviance df.resid 
##    387.1    403.4   -187.5    375.1      107 
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -1.55254 -0.39337  0.03275  0.37793  1.69016 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  subj_id  (Intercept) 1.1179   1.0573  
##  Residual             0.6105   0.7813  
## Number of obs: 113, groups:  subj_id, 102
## 
## Fixed effects:
##                  Estimate Std. Error       df t value Pr(>|t|)    
## (Intercept)       13.6886     0.6657 106.4622  20.562   <2e-16 ***
## mri_age_y          0.2954     0.1616 107.8036   1.827   0.0704 .  
## female             0.2522     0.9254 110.1751   0.273   0.7857    
## mri_age_y:female  -0.0757     0.2178 109.0452  -0.348   0.7288    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) mr_g_y female
## mri_age_y   -0.964              
## female      -0.719  0.694       
## mri_g_y:fml  0.716 -0.742 -0.960
anova(acg_naa_age,acg_naa_ageXsex)
## Data: acg_dat
## Models:
## acg_naa_age: NAA_conc_molal ~ mri_age_y + (1 | subj_id)
## acg_naa_ageXsex: NAA_conc_molal ~ mri_age_y + female + mri_age_y:female + (1 | subj_id)
##                 npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)
## acg_naa_age        4 383.25 394.16 -187.62   375.25                     
## acg_naa_ageXsex    6 387.08 403.45 -187.54   375.08 0.1677  2     0.9196
#including sex in model does not significantly improve fit 

Cre

acg_cre_age<- lmer(Cre_PCr_conc_molal ~ mri_age_y + (1 | subj_id), data=acg_dat, REML=FALSE)
summary(acg_cre_age)
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
##   method [lmerModLmerTest]
## Formula: Cre_PCr_conc_molal ~ mri_age_y + (1 | subj_id)
##    Data: acg_dat
## 
##      AIC      BIC   logLik deviance df.resid 
##    340.9    351.8   -166.5    332.9      109 
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.12075 -0.39494  0.02556  0.37255  1.62021 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  subj_id  (Intercept) 0.6867   0.8287  
##  Residual             0.4877   0.6983  
## Number of obs: 113, groups:  subj_id, 102
## 
## Fixed effects:
##              Estimate Std. Error        df t value Pr(>|t|)    
## (Intercept)  11.67366    0.38192 110.29328  30.566   <2e-16 ***
## mri_age_y     0.09221    0.08992 107.78610   1.026    0.307    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##           (Intr)
## mri_age_y -0.960
acg_cre_age_sex<- lmer(Cre_PCr_conc_molal ~ mri_age_y + female + (1 | subj_id), data=acg_dat,  REML=FALSE)
summary(acg_cre_age_sex)
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
##   method [lmerModLmerTest]
## Formula: Cre_PCr_conc_molal ~ mri_age_y + female + (1 | subj_id)
##    Data: acg_dat
## 
##      AIC      BIC   logLik deviance df.resid 
##    342.9    356.5   -166.4    332.9      108 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.1065 -0.3932  0.0080  0.3605  1.6020 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  subj_id  (Intercept) 0.6866   0.8286  
##  Residual             0.4871   0.6979  
## Number of obs: 113, groups:  subj_id, 102
## 
## Fixed effects:
##              Estimate Std. Error        df t value Pr(>|t|)    
## (Intercept)  11.65746    0.38678 111.48623  30.139   <2e-16 ***
## mri_age_y     0.08988    0.09034 107.34558   0.995    0.322    
## female        0.05586    0.21445  96.67410   0.260    0.795    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##           (Intr) mr_g_y
## mri_age_y -0.927       
## female    -0.160 -0.100
anova(acg_cre_age,acg_cre_age_sex)
## Data: acg_dat
## Models:
## acg_cre_age: Cre_PCr_conc_molal ~ mri_age_y + (1 | subj_id)
## acg_cre_age_sex: Cre_PCr_conc_molal ~ mri_age_y + female + (1 | subj_id)
##                 npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)
## acg_cre_age        4 340.94 351.85 -166.47   332.94                     
## acg_cre_age_sex    5 342.87 356.51 -166.44   332.87 0.0678  1     0.7945
acg_cre_ageXsex<- lmer(Cre_PCr_conc_molal ~ mri_age_y + female + mri_age_y:female + (1 | subj_id), data=acg_dat, REML=FALSE)
summary(acg_cre_ageXsex)
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
##   method [lmerModLmerTest]
## Formula: Cre_PCr_conc_molal ~ mri_age_y + female + mri_age_y:female +  
##     (1 | subj_id)
##    Data: acg_dat
## 
##      AIC      BIC   logLik deviance df.resid 
##    343.6    360.0   -165.8    331.6      107 
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.11740 -0.37592  0.04868  0.38312  1.59263 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  subj_id  (Intercept) 0.7009   0.8372  
##  Residual             0.4638   0.6810  
## Number of obs: 113, groups:  subj_id, 102
## 
## Fixed effects:
##                   Estimate Std. Error        df t value Pr(>|t|)    
## (Intercept)       12.09093    0.54699 105.10134  22.104   <2e-16 ***
## mri_age_y         -0.01931    0.13288 106.80701  -0.145    0.885    
## female            -0.77294    0.76525 111.06817  -1.010    0.315    
## mri_age_y:female   0.20325    0.18016 110.02199   1.128    0.262    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) mr_g_y female
## mri_age_y   -0.965              
## female      -0.715  0.689       
## mri_g_y:fml  0.711 -0.738 -0.960
anova(acg_cre_age,acg_cre_ageXsex)
## Data: acg_dat
## Models:
## acg_cre_age: Cre_PCr_conc_molal ~ mri_age_y + (1 | subj_id)
## acg_cre_ageXsex: Cre_PCr_conc_molal ~ mri_age_y + female + mri_age_y:female + (1 | subj_id)
##                 npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)
## acg_cre_age        4 340.94 351.85 -166.47   332.94                     
## acg_cre_ageXsex    6 343.61 359.98 -165.81   331.61 1.3273  2      0.515

Cho

acg_cho_age<- lmer(Cho_GPC_PCh_conc_molal ~ mri_age_y + (1 | subj_id), data=acg_dat,  REML=FALSE)
summary(acg_cho_age)
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
##   method [lmerModLmerTest]
## Formula: Cho_GPC_PCh_conc_molal ~ mri_age_y + (1 | subj_id)
##    Data: acg_dat
## 
##      AIC      BIC   logLik deviance df.resid 
##     78.1     89.0    -35.0     70.1      108 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.1325 -0.5979 -0.1286  0.6106  2.4400 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  subj_id  (Intercept) 0.03204  0.1790  
##  Residual             0.07887  0.2808  
## Number of obs: 112, groups:  subj_id, 101
## 
## Fixed effects:
##              Estimate Std. Error        df t value Pr(>|t|)    
## (Intercept)   2.74592    0.12031 111.92574   22.82   <2e-16 ***
## mri_age_y    -0.06423    0.02842 111.94035   -2.26   0.0258 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##           (Intr)
## mri_age_y -0.963
acg_cho_age_sex<- lmer(Cho_GPC_PCh_conc_molal ~ mri_age_y + female + (1 | subj_id), data=acg_dat,  REML=FALSE)
summary(acg_cho_age_sex)
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
##   method [lmerModLmerTest]
## Formula: Cho_GPC_PCh_conc_molal ~ mri_age_y + female + (1 | subj_id)
##    Data: acg_dat
## 
##      AIC      BIC   logLik deviance df.resid 
##     78.2     91.8    -34.1     68.2      107 
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.17836 -0.63945 -0.09041  0.54644  2.50032 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  subj_id  (Intercept) 0.03214  0.1793  
##  Residual             0.07694  0.2774  
## Number of obs: 112, groups:  subj_id, 101
## 
## Fixed effects:
##              Estimate Std. Error        df t value Pr(>|t|)    
## (Intercept)   2.77074    0.12063 111.75493  22.970   <2e-16 ***
## mri_age_y    -0.06004    0.02834 111.86929  -2.118   0.0364 *  
## female       -0.09017    0.06493 103.29750  -1.389   0.1679    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##           (Intr) mr_g_y
## mri_age_y -0.931       
## female    -0.148 -0.107
anova(acg_cho_age,acg_cho_age_sex)
## Data: acg_dat
## Models:
## acg_cho_age: Cho_GPC_PCh_conc_molal ~ mri_age_y + (1 | subj_id)
## acg_cho_age_sex: Cho_GPC_PCh_conc_molal ~ mri_age_y + female + (1 | subj_id)
##                 npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)
## acg_cho_age        4 78.077 88.951 -35.039   70.077                     
## acg_cho_age_sex    5 78.166 91.759 -34.083   68.166 1.9111  1     0.1668
acg_cho_ageXsex<- lmer(Cho_GPC_PCh_conc_molal ~ mri_age_y + female + mri_age_y:female + (1 | subj_id), data=acg_dat,  REML=FALSE)
summary(acg_cho_ageXsex)
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
##   method [lmerModLmerTest]
## Formula: Cho_GPC_PCh_conc_molal ~ mri_age_y + female + mri_age_y:female +  
##     (1 | subj_id)
##    Data: acg_dat
## 
##      AIC      BIC   logLik deviance df.resid 
##     79.7     96.1    -33.9     67.7      106 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.1308 -0.6109 -0.1102  0.5606  2.5638 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  subj_id  (Intercept) 0.03297  0.1816  
##  Residual             0.07579  0.2753  
## Number of obs: 112, groups:  subj_id, 101
## 
## Fixed effects:
##                   Estimate Std. Error        df t value Pr(>|t|)    
## (Intercept)        2.69546    0.16726 107.68134  16.116   <2e-16 ***
## mri_age_y         -0.04107    0.04072 109.53005  -1.008    0.315    
## female             0.05961    0.23991 111.98464   0.248    0.804    
## mri_age_y:female  -0.03674    0.05661 111.93901  -0.649    0.518    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) mr_g_y female
## mri_age_y   -0.965              
## female      -0.697  0.673       
## mri_g_y:fml  0.694 -0.719 -0.963
anova(acg_cho_age,acg_cho_ageXsex)
## Data: acg_dat
## Models:
## acg_cho_age: Cho_GPC_PCh_conc_molal ~ mri_age_y + (1 | subj_id)
## acg_cho_ageXsex: Cho_GPC_PCh_conc_molal ~ mri_age_y + female + mri_age_y:female + (1 | subj_id)
##                 npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)
## acg_cho_age        4 78.077 88.951 -35.039   70.077                     
## acg_cho_ageXsex    6 79.747 96.058 -33.874   67.747 2.3297  2      0.312
#including sex in model does not significantly improve fit 

Ins

acg_ins_age<- lmer(Ins_conc_molal ~ mri_age_y + (1 | subj_id), data=acg_dat, REML=FALSE)
summary(acg_ins_age)
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
##   method [lmerModLmerTest]
## Formula: Ins_conc_molal ~ mri_age_y + (1 | subj_id)
##    Data: acg_dat
## 
##      AIC      BIC   logLik deviance df.resid 
##    416.5    427.3   -204.2    408.5      108 
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.51194 -0.47558  0.09217  0.51196  1.57868 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  subj_id  (Intercept) 0.9268   0.9627  
##  Residual             1.3749   1.1726  
## Number of obs: 112, groups:  subj_id, 101
## 
## Fixed effects:
##             Estimate Std. Error       df t value Pr(>|t|)    
## (Intercept)   8.4742     0.5464 111.9652  15.509   <2e-16 ***
## mri_age_y    -0.2049     0.1294 111.4483  -1.583    0.116    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##           (Intr)
## mri_age_y -0.962
acg_ins_age_sex<- lmer(Ins_conc_molal ~ mri_age_y + female + (1 | subj_id), data=acg_dat, REML=FALSE)
summary(acg_ins_age_sex)
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
##   method [lmerModLmerTest]
## Formula: Ins_conc_molal ~ mri_age_y + female + (1 | subj_id)
##    Data: acg_dat
## 
##      AIC      BIC   logLik deviance df.resid 
##    417.1    430.7   -203.5    407.1      107 
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.47103 -0.48748  0.05729  0.54401  1.51171 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  subj_id  (Intercept) 0.8891   0.9429  
##  Residual             1.3819   1.1755  
## Number of obs: 112, groups:  subj_id, 101
## 
## Fixed effects:
##             Estimate Std. Error       df t value Pr(>|t|)    
## (Intercept)   8.5595     0.5485 111.9724  15.605   <2e-16 ***
## mri_age_y    -0.1859     0.1296 111.3196  -1.435    0.154    
## female       -0.3489     0.2981  99.5526  -1.171    0.245    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##           (Intr) mr_g_y
## mri_age_y -0.930       
## female    -0.140 -0.118
anova(acg_ins_age,acg_ins_age_sex)
## Data: acg_dat
## Models:
## acg_ins_age: Ins_conc_molal ~ mri_age_y + (1 | subj_id)
## acg_ins_age_sex: Ins_conc_molal ~ mri_age_y + female + (1 | subj_id)
##                 npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)
## acg_ins_age        4 416.45 427.33 -204.23   408.45                     
## acg_ins_age_sex    5 417.09 430.69 -203.55   407.09 1.3598  1     0.2436
acg_ins_ageXsex<- lmer(Ins_conc_molal ~ mri_age_y + female + mri_age_y:female + (1 | subj_id), data=acg_dat, REML=FALSE)
summary(acg_ins_ageXsex)
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
##   method [lmerModLmerTest]
## Formula: Ins_conc_molal ~ mri_age_y + female + mri_age_y:female + (1 |  
##     subj_id)
##    Data: acg_dat
## 
##      AIC      BIC   logLik deviance df.resid 
##    418.2    434.5   -203.1    406.2      106 
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.55603 -0.53421  0.07378  0.57588  1.65390 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  subj_id  (Intercept) 0.7615   0.8727  
##  Residual             1.4790   1.2162  
## Number of obs: 112, groups:  subj_id, 101
## 
## Fixed effects:
##                  Estimate Std. Error       df t value Pr(>|t|)    
## (Intercept)        9.0615     0.7667 105.7180  11.819   <2e-16 ***
## mri_age_y         -0.3138     0.1882 108.1469  -1.667   0.0983 .  
## female            -1.3693     1.0920 111.9981  -1.254   0.2125    
## mri_age_y:female   0.2511     0.2588 111.8723   0.970   0.3340    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) mr_g_y female
## mri_age_y   -0.965              
## female      -0.702  0.678       
## mri_g_y:fml  0.702 -0.727 -0.963
anova(acg_ins_age,acg_ins_ageXsex)
## Data: acg_dat
## Models:
## acg_ins_age: Ins_conc_molal ~ mri_age_y + (1 | subj_id)
## acg_ins_ageXsex: Ins_conc_molal ~ mri_age_y + female + mri_age_y:female + (1 | subj_id)
##                 npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)
## acg_ins_age        4 416.45 427.33 -204.23   408.45                     
## acg_ins_ageXsex    6 418.20 434.51 -203.10   406.20 2.2549  2     0.3239

Glx

acg_glx_age<- lmer(Glu_Gln_conc_molal ~ mri_age_y + (1 | subj_id), data=acg_dat, REML=FALSE)
summary(acg_glx_age)
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
##   method [lmerModLmerTest]
## Formula: Glu_Gln_conc_molal ~ mri_age_y + (1 | subj_id)
##    Data: acg_dat
## 
##      AIC      BIC   logLik deviance df.resid 
##    541.1    552.0   -266.5    533.1      109 
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -1.48535 -0.37243 -0.00651  0.43846  1.68169 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  subj_id  (Intercept) 4.706    2.169   
##  Residual             2.330    1.526   
## Number of obs: 113, groups:  subj_id, 102
## 
## Fixed effects:
##             Estimate Std. Error       df t value Pr(>|t|)    
## (Intercept)  25.7502     0.9206 106.5717  27.972   <2e-16 ***
## mri_age_y     0.2612     0.2163 103.0785   1.207     0.23    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##           (Intr)
## mri_age_y -0.959
acg_glx_age_sex<- lmer(Glu_Gln_conc_molal ~ mri_age_y + female + (1 | subj_id), data=acg_dat, REML=FALSE)
summary(acg_glx_age_sex)
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
##   method [lmerModLmerTest]
## Formula: Glu_Gln_conc_molal ~ mri_age_y + female + (1 | subj_id)
##    Data: acg_dat
## 
##      AIC      BIC   logLik deviance df.resid 
##    543.0    556.6   -266.5    533.0      108 
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -1.47879 -0.39442 -0.02857  0.42851  1.70104 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  subj_id  (Intercept) 4.692    2.166   
##  Residual             2.335    1.528   
## Number of obs: 113, groups:  subj_id, 102
## 
## Fixed effects:
##             Estimate Std. Error       df t value Pr(>|t|)    
## (Intercept)  25.7002     0.9334 108.5734  27.535   <2e-16 ***
## mri_age_y     0.2542     0.2173 102.5984   1.170    0.245    
## female        0.1700     0.5257  98.8828   0.323    0.747    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##           (Intr) mr_g_y
## mri_age_y -0.925       
## female    -0.167 -0.097
anova(acg_glx_age,acg_glx_age_sex)
## Data: acg_dat
## Models:
## acg_glx_age: Glu_Gln_conc_molal ~ mri_age_y + (1 | subj_id)
## acg_glx_age_sex: Glu_Gln_conc_molal ~ mri_age_y + female + (1 | subj_id)
##                 npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)
## acg_glx_age        4 541.09 552.00 -266.55   533.09                     
## acg_glx_age_sex    5 542.99 556.62 -266.49   532.99 0.1045  1     0.7465
acg_glx_ageXsex<- lmer(Glu_Gln_conc_molal ~ mri_age_y + female + mri_age_y:female + (1 | subj_id), data=acg_dat, REML=FALSE)
summary(acg_glx_ageXsex)
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
##   method [lmerModLmerTest]
## Formula: Glu_Gln_conc_molal ~ mri_age_y + female + mri_age_y:female +  
##     (1 | subj_id)
##    Data: acg_dat
## 
##      AIC      BIC   logLik deviance df.resid 
##    542.5    558.9   -265.2    530.5      107 
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -1.53333 -0.43763  0.07324  0.38166  1.63409 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  subj_id  (Intercept) 4.573    2.139   
##  Residual             2.297    1.515   
## Number of obs: 113, groups:  subj_id, 102
## 
## Fixed effects:
##                  Estimate Std. Error       df t value Pr(>|t|)    
## (Intercept)       27.2150     1.3265 106.5825  20.516   <2e-16 ***
## mri_age_y         -0.1271     0.3220 107.8777  -0.395    0.694    
## female            -2.6328     1.8384 109.4802  -1.432    0.155    
## mri_age_y:female   0.6875     0.4325 108.2254   1.590    0.115    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) mr_g_y female
## mri_age_y   -0.964              
## female      -0.722  0.696       
## mri_g_y:fml  0.718 -0.745 -0.959
anova(acg_glx_age,acg_glx_ageXsex)
## Data: acg_dat
## Models:
## acg_glx_age: Glu_Gln_conc_molal ~ mri_age_y + (1 | subj_id)
## acg_glx_ageXsex: Glu_Gln_conc_molal ~ mri_age_y + female + mri_age_y:female + (1 | subj_id)
##                 npar    AIC    BIC  logLik deviance Chisq Df Pr(>Chisq)
## acg_glx_age        4 541.09 552.00 -266.55   533.09                    
## acg_glx_ageXsex    6 542.49 558.85 -265.25   530.49 2.603  2     0.2721

Plot results ACG

Create scatter plot with regression line based on model predictions NAA

#pull intercept and slope from model for regression line
fixef.acg_naa<-fixef(acg_naa_age)

#plot
plot_acg_naa_age<-ggplot(acg_dat, aes(x=mri_age_y, y=NAA_conc_molal,color=as.factor(female)))
plot_acg_naa_age + 
  geom_point() +
  geom_line(aes(group = as.factor(subj_id))) +
  labs(title = "Anterior Cingulate", x = "Age (years)", y = "tNAA (molal)") +
  scale_color_discrete(name = "Sex", labels = c("Male", "Female")) +
  geom_abline(aes(intercept = fixef.acg_naa[[1]], slope= fixef.acg_naa[[2]])) + 
  theme_classic()

ggsave("/Users/meaghan/OneDrive - University of Calgary/Preschool_data/MRS/results/acg_naa_age.png")
## Saving 7 x 5 in image

Cho

#pull intercept and slope from model for regression line
fixef.acg_cho<-fixef(acg_cho_age)

#plot
plot_acg_cho_age<-ggplot(acg_dat, aes(x=mri_age_y, y=Cho_GPC_PCh_conc_molal,color=as.factor(female)))
plot_acg_cho_age + 
  geom_point() +
  geom_line(aes(group = as.factor(subj_id))) +
  labs(title = "Anterior Cingulate", x = "Age (years)", y = "tCho (molal)") +
  scale_color_discrete(name = "Sex", labels = c("Male", "Female")) +
  geom_abline(aes(intercept = fixef.acg_cho[[1]], slope= fixef.acg_cho[[2]])) + 
  theme_classic()
## Warning: Removed 1 rows containing missing values (geom_point).
## Warning: Removed 1 row(s) containing missing values (geom_path).

ggsave("/Users/meaghan/OneDrive - University of Calgary/Preschool_data/MRS/results/acg_cho_age.png")
## Saving 7 x 5 in image
## Warning: Removed 1 rows containing missing values (geom_point).

## Warning: Removed 1 row(s) containing missing values (geom_path).

Ins

#pull intercept and slope from model for regression line
fixef.acg_ins<-fixef(acg_ins_age)

#plot
plot_acg_ins_age<-ggplot(acg_dat, aes(x=mri_age_y, y=Ins_conc_molal,color=as.factor(female)))
plot_acg_ins_age + 
  geom_point() +
  geom_line(aes(group = as.factor(subj_id))) +
  labs(title = "Anterior Cingulate", x = "Age (years)", y = "Ins (molal)") +
  scale_color_discrete(name = "Sex", labels = c("Male", "Female")) +
  geom_abline(aes(intercept = fixef.acg_ins[[1]], slope= fixef.acg_ins[[2]])) + 
  theme_classic()
## Warning: Removed 1 rows containing missing values (geom_point).
## Warning: Removed 1 row(s) containing missing values (geom_path).

ggsave("/Users/meaghan/OneDrive - University of Calgary/Preschool_data/MRS/results/acg_ins_age.png")
## Saving 7 x 5 in image
## Warning: Removed 1 rows containing missing values (geom_point).

## Warning: Removed 1 row(s) containing missing values (geom_path).

LME modelling LAG

run mixed-effects models for lag with participant as random effect, age, sex and age x sex interaction as fixed effects NAA, substantial longitudinal data so we include random slope run on lag_data_molal_outliers_removed_n320.csv

#basic model effect of age on intercept
lag_naa_age<- lmer(NAA_conc_molal ~ mri_age_y + (1 | subj_id), data=lag_dat, REML=FALSE)
summary(lag_naa_age)
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
##   method [lmerModLmerTest]
## Formula: NAA_conc_molal ~ mri_age_y + (1 | subj_id)
##    Data: lag_dat
## 
##      AIC      BIC   logLik deviance df.resid 
##      896      911     -444      888      314 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.4313 -0.5792 -0.0185  0.5440  3.6263 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  subj_id  (Intercept) 0.1641   0.4051  
##  Residual             0.8272   0.9095  
## Number of obs: 318, groups:  subj_id, 95
## 
## Fixed effects:
##              Estimate Std. Error        df t value Pr(>|t|)    
## (Intercept)  13.82375    0.18588 309.72511  74.369  < 2e-16 ***
## mri_age_y     0.09699    0.02902 316.51634   3.343  0.00093 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##           (Intr)
## mri_age_y -0.929
lag_naa_age_sex<- lmer(NAA_conc_molal ~ mri_age_y + female + (1 | subj_id), data=lag_dat, REML=FALSE)
summary(lag_naa_age_sex)
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
##   method [lmerModLmerTest]
## Formula: NAA_conc_molal ~ mri_age_y + female + (1 | subj_id)
##    Data: lag_dat
## 
##      AIC      BIC   logLik deviance df.resid 
##    898.0    916.8   -444.0    888.0      313 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.4286 -0.5791 -0.0185  0.5452  3.6231 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  subj_id  (Intercept) 0.1637   0.4046  
##  Residual             0.8275   0.9097  
## Number of obs: 318, groups:  subj_id, 95
## 
## Fixed effects:
##               Estimate Std. Error         df t value Pr(>|t|)    
## (Intercept)  13.827819   0.194854 290.918032  70.965  < 2e-16 ***
## mri_age_y     0.097086   0.029038 316.589446   3.343 0.000927 ***
## female       -0.009878   0.137557  93.706753  -0.072 0.942903    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##           (Intr) mr_g_y
## mri_age_y -0.874       
## female    -0.300 -0.039
anova(lag_naa_age,lag_naa_age_sex)
## Data: lag_dat
## Models:
## lag_naa_age: NAA_conc_molal ~ mri_age_y + (1 | subj_id)
## lag_naa_age_sex: NAA_conc_molal ~ mri_age_y + female + (1 | subj_id)
##                 npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)
## lag_naa_age        4 895.96 911.01 -443.98   887.96                     
## lag_naa_age_sex    5 897.95 916.76 -443.98   887.95 0.0051  1      0.943
lag_naa_ageXsex<- lmer(NAA_conc_molal ~ mri_age_y + female + mri_age_y:female + (1 | subj_id), data=lag_dat, REML=FALSE)
summary(lag_naa_ageXsex)
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
##   method [lmerModLmerTest]
## Formula: NAA_conc_molal ~ mri_age_y + female + mri_age_y:female + (1 |  
##     subj_id)
##    Data: lag_dat
## 
##      AIC      BIC   logLik deviance df.resid 
##    899.9    922.5   -444.0    887.9      312 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.4372 -0.5859 -0.0225  0.5510  3.6185 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  subj_id  (Intercept) 0.1645   0.4056  
##  Residual             0.8268   0.9093  
## Number of obs: 318, groups:  subj_id, 95
## 
## Fixed effects:
##                   Estimate Std. Error        df t value Pr(>|t|)    
## (Intercept)       13.79473    0.24812 314.11810  55.598  < 2e-16 ***
## mri_age_y          0.10272    0.03910 313.83715   2.627  0.00903 ** 
## female             0.06577    0.37459 308.18552   0.176  0.86075    
## mri_age_y:female  -0.01265    0.05838 316.98288  -0.217  0.82864    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) mr_g_y female
## mri_age_y   -0.924              
## female      -0.662  0.612       
## mri_g_y:fml  0.619 -0.670 -0.930
anova(lag_naa_age,lag_naa_ageXsex)
## Data: lag_dat
## Models:
## lag_naa_age: NAA_conc_molal ~ mri_age_y + (1 | subj_id)
## lag_naa_ageXsex: NAA_conc_molal ~ mri_age_y + female + mri_age_y:female + (1 | subj_id)
##                 npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)
## lag_naa_age        4 895.96 911.01 -443.98   887.96                     
## lag_naa_ageXsex    6 899.91 922.48 -443.95   887.91 0.0518  2     0.9744
#including sex in model does not significantly improve fit 

#test quadratic effect of age
#create age^2  variable
lag_dat$mri_age_y2 <- (lag_dat$mri_age_y)^2

lag_naa_age2<- lmer(NAA_conc_molal ~ mri_age_y + mri_age_y2 + (1 | subj_id), data=lag_dat, REML=FALSE)
summary(lag_naa_age2) #age^2 p=.0447
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
##   method [lmerModLmerTest]
## Formula: NAA_conc_molal ~ mri_age_y + mri_age_y2 + (1 | subj_id)
##    Data: lag_dat
## 
##      AIC      BIC   logLik deviance df.resid 
##    893.9    912.7   -442.0    883.9      313 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.5451 -0.5819 -0.0270  0.4941  3.5330 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  subj_id  (Intercept) 0.1700   0.4123  
##  Residual             0.8119   0.9011  
## Number of obs: 318, groups:  subj_id, 95
## 
## Fixed effects:
##              Estimate Std. Error        df t value Pr(>|t|)    
## (Intercept)  12.70181    0.58588 310.95949  21.680   <2e-16 ***
## mri_age_y     0.47130    0.18749 303.10334   2.514   0.0125 *  
## mri_age_y2   -0.02836    0.01402 297.51463  -2.023   0.0440 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##            (Intr) mr_g_y
## mri_age_y  -0.983       
## mri_age_y2  0.949 -0.988
anova(lag_naa_age, lag_naa_age2) #fit improved p=0.44
## Data: lag_dat
## Models:
## lag_naa_age: NAA_conc_molal ~ mri_age_y + (1 | subj_id)
## lag_naa_age2: NAA_conc_molal ~ mri_age_y + mri_age_y2 + (1 | subj_id)
##              npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)  
## lag_naa_age     4 895.96 911.01 -443.98   887.96                       
## lag_naa_age2    5 893.91 912.72 -441.96   883.91 4.0455  1    0.04429 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Cre

lag_cre_age<- lmer(Cre_PCr_conc_molal ~ mri_age_y + (1 | subj_id), data=lag_dat, REML=FALSE)
summary(lag_cre_age)
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
##   method [lmerModLmerTest]
## Formula: Cre_PCr_conc_molal ~ mri_age_y + (1 | subj_id)
##    Data: lag_dat
## 
##      AIC      BIC   logLik deviance df.resid 
##    848.2    863.3   -420.1    840.2      316 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.9764 -0.5307 -0.1261  0.5525  4.1548 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  subj_id  (Intercept) 0.08477  0.2912  
##  Residual             0.73683  0.8584  
## Number of obs: 320, groups:  subj_id, 95
## 
## Fixed effects:
##              Estimate Std. Error        df t value Pr(>|t|)    
## (Intercept)  10.49912    0.16805 308.49608  62.476   <2e-16 ***
## mri_age_y    -0.04183    0.02646 319.87112  -1.581    0.115    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##           (Intr)
## mri_age_y -0.937
lag_cre_age_sex<- lmer(Cre_PCr_conc_molal ~ mri_age_y + female + (1 | subj_id), data=lag_dat, REML=FALSE)
summary(lag_cre_age_sex)
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
##   method [lmerModLmerTest]
## Formula: Cre_PCr_conc_molal ~ mri_age_y + female + (1 | subj_id)
##    Data: lag_dat
## 
##      AIC      BIC   logLik deviance df.resid 
##    849.7    868.6   -419.9    839.7      315 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.9528 -0.5658 -0.1208  0.5434  4.1925 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  subj_id  (Intercept) 0.08478  0.2912  
##  Residual             0.73557  0.8577  
## Number of obs: 320, groups:  subj_id, 95
## 
## Fixed effects:
##              Estimate Std. Error        df t value Pr(>|t|)    
## (Intercept)  10.46482    0.17478 289.38001  59.873   <2e-16 ***
## mri_age_y    -0.04267    0.02646 319.88902  -1.612    0.108    
## female        0.08318    0.11739  84.73092   0.709    0.481    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##           (Intr) mr_g_y
## mri_age_y -0.887       
## female    -0.277 -0.044
anova(lag_cre_age,lag_cre_age_sex)
## Data: lag_dat
## Models:
## lag_cre_age: Cre_PCr_conc_molal ~ mri_age_y + (1 | subj_id)
## lag_cre_age_sex: Cre_PCr_conc_molal ~ mri_age_y + female + (1 | subj_id)
##                 npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)
## lag_cre_age        4 848.23 863.31 -420.12   840.23                     
## lag_cre_age_sex    5 849.73 868.57 -419.87   839.73 0.5017  1     0.4788
lag_cre_ageXsex<- lmer(Cre_PCr_conc_molal ~ mri_age_y + female + mri_age_y:female + (1 | subj_id), data=lag_dat, REML=FALSE)
summary(lag_cre_ageXsex)
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
##   method [lmerModLmerTest]
## Formula: Cre_PCr_conc_molal ~ mri_age_y + female + mri_age_y:female +  
##     (1 | subj_id)
##    Data: lag_dat
## 
##      AIC      BIC   logLik deviance df.resid 
##    850.0    872.6   -419.0    838.0      314 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.9579 -0.5653 -0.0910  0.5718  4.2323 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  subj_id  (Intercept) 0.08177  0.2860  
##  Residual             0.73339  0.8564  
## Number of obs: 320, groups:  subj_id, 95
## 
## Fixed effects:
##                   Estimate Std. Error        df t value Pr(>|t|)    
## (Intercept)       10.65304    0.22484 314.11103  47.381   <2e-16 ***
## mri_age_y         -0.07476    0.03585 318.90184  -2.086   0.0378 *  
## female            -0.33761    0.33688 307.33953  -1.002   0.3171    
## mri_age_y:female   0.07047    0.05297 319.97938   1.330   0.1844    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) mr_g_y female
## mri_age_y   -0.935              
## female      -0.667  0.624       
## mri_g_y:fml  0.632 -0.677 -0.938
anova(lag_cre_age,lag_cre_ageXsex)
## Data: lag_dat
## Models:
## lag_cre_age: Cre_PCr_conc_molal ~ mri_age_y + (1 | subj_id)
## lag_cre_ageXsex: Cre_PCr_conc_molal ~ mri_age_y + female + mri_age_y:female + (1 | subj_id)
##                 npar    AIC    BIC  logLik deviance Chisq Df Pr(>Chisq)
## lag_cre_age        4 848.23 863.31 -420.12   840.23                    
## lag_cre_ageXsex    6 849.97 872.58 -418.98   837.97 2.263  2     0.3226

Cho

lag_cho_age<- lmer(Cho_GPC_PCh_conc_molal ~ mri_age_y + (1 | subj_id), data=lag_dat, REML=FALSE)
summary(lag_cho_age)
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
##   method [lmerModLmerTest]
## Formula: Cho_GPC_PCh_conc_molal ~ mri_age_y + (1 | subj_id)
##    Data: lag_dat
## 
##      AIC      BIC   logLik deviance df.resid 
##     51.3     66.4    -21.7     43.3      316 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.6539 -0.5810 -0.0562  0.4659  4.1202 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  subj_id  (Intercept) 0.0149   0.1221  
##  Residual             0.0561   0.2368  
## Number of obs: 320, groups:  subj_id, 95
## 
## Fixed effects:
##               Estimate Std. Error         df t value Pr(>|t|)    
## (Intercept)   2.346670   0.049038 310.144986  47.854  < 2e-16 ***
## mri_age_y    -0.023315   0.007581 315.773363  -3.075  0.00229 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##           (Intr)
## mri_age_y -0.921
lag_cho_age_sex<- lmer(Cho_GPC_PCh_conc_molal ~ mri_age_y + female + (1 | subj_id), data=lag_dat, REML=FALSE)
summary(lag_cho_age_sex)
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
##   method [lmerModLmerTest]
## Formula: Cho_GPC_PCh_conc_molal ~ mri_age_y + female + (1 | subj_id)
##    Data: lag_dat
## 
##      AIC      BIC   logLik deviance df.resid 
##     51.6     70.5    -20.8     41.6      315 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.6317 -0.5893 -0.0587  0.4688  4.1717 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  subj_id  (Intercept) 0.01469  0.1212  
##  Residual             0.05586  0.2364  
## Number of obs: 320, groups:  subj_id, 95
## 
## Fixed effects:
##              Estimate Std. Error        df t value Pr(>|t|)    
## (Intercept)   2.36786    0.05147 283.09492  46.006  < 2e-16 ***
## mri_age_y    -0.02286    0.00757 315.96555  -3.019  0.00274 ** 
## female       -0.05024    0.03808  80.73098  -1.319  0.19083    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##           (Intr) mr_g_y
## mri_age_y -0.860       
## female    -0.312 -0.046
anova(lag_cho_age,lag_cho_age_sex)
## Data: lag_dat
## Models:
## lag_cho_age: Cho_GPC_PCh_conc_molal ~ mri_age_y + (1 | subj_id)
## lag_cho_age_sex: Cho_GPC_PCh_conc_molal ~ mri_age_y + female + (1 | subj_id)
##                 npar    AIC    BIC  logLik deviance Chisq Df Pr(>Chisq)
## lag_cho_age        4 51.345 66.419 -21.673   43.345                    
## lag_cho_age_sex    5 51.610 70.452 -20.805   41.610 1.735  1     0.1878
lag_cho_ageXsex<- lmer(Cho_GPC_PCh_conc_molal ~ mri_age_y + female + mri_age_y:female + (1 | subj_id), data=lag_dat, REML=FALSE)
summary(lag_cho_ageXsex)
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
##   method [lmerModLmerTest]
## Formula: Cho_GPC_PCh_conc_molal ~ mri_age_y + female + mri_age_y:female +  
##     (1 | subj_id)
##    Data: lag_dat
## 
##      AIC      BIC   logLik deviance df.resid 
##     49.9     72.5    -18.9     37.9      314 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.6639 -0.5411 -0.0413  0.4642  4.0526 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  subj_id  (Intercept) 0.01452  0.1205  
##  Residual             0.05521  0.2350  
## Number of obs: 320, groups:  subj_id, 95
## 
## Fixed effects:
##                   Estimate Std. Error        df t value Pr(>|t|)    
## (Intercept)        2.28998    0.06505 314.68723  35.206   <2e-16 ***
## mri_age_y         -0.00954    0.01019 312.12139  -0.936   0.3497    
## female             0.12488    0.09791 309.31059   1.275   0.2031    
## mri_age_y:female  -0.02931    0.01511 316.65469  -1.939   0.0533 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) mr_g_y female
## mri_age_y   -0.916              
## female      -0.664  0.609       
## mri_g_y:fml  0.617 -0.674 -0.922
anova(lag_cho_age,lag_cho_ageXsex)
## Data: lag_dat
## Models:
## lag_cho_age: Cho_GPC_PCh_conc_molal ~ mri_age_y + (1 | subj_id)
## lag_cho_ageXsex: Cho_GPC_PCh_conc_molal ~ mri_age_y + female + mri_age_y:female + (1 | subj_id)
##                 npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)  
## lag_cho_age        4 51.345 66.419 -21.673   43.345                       
## lag_cho_ageXsex    6 49.871 72.481 -18.936   37.871 5.4741  2    0.06476 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#including sex in model does not significantly improve fit 

Ins

lag_ins_age<- lmer(Ins_conc_molal ~ mri_age_y + (1 | subj_id), data=lag_dat, REML=FALSE)
summary(lag_ins_age)
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
##   method [lmerModLmerTest]
## Formula: Ins_conc_molal ~ mri_age_y + (1 | subj_id)
##    Data: lag_dat
## 
##      AIC      BIC   logLik deviance df.resid 
##   1013.3   1028.3   -502.6   1005.3      315 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.3984 -0.5594 -0.0111  0.5911  4.6516 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  subj_id  (Intercept) 0.1768   0.4205  
##  Residual             1.2227   1.1057  
## Number of obs: 319, groups:  subj_id, 95
## 
## Fixed effects:
##              Estimate Std. Error        df t value Pr(>|t|)    
## (Intercept)   6.28473    0.21938 307.28604  28.648   <2e-16 ***
## mri_age_y     0.01632    0.03444 318.45543   0.474    0.636    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##           (Intr)
## mri_age_y -0.934
lag_ins_age_sex<- lmer(Ins_conc_molal ~ mri_age_y + female + (1 | subj_id), data=lag_dat, REML=FALSE)
summary(lag_ins_age_sex)
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
##   method [lmerModLmerTest]
## Formula: Ins_conc_molal ~ mri_age_y + female + (1 | subj_id)
##    Data: lag_dat
## 
##      AIC      BIC   logLik deviance df.resid 
##   1015.3   1034.1   -502.6   1005.3      314 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.4003 -0.5588 -0.0120  0.5901  4.6506 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  subj_id  (Intercept) 0.1765   0.4201  
##  Residual             1.2229   1.1058  
## Number of obs: 319, groups:  subj_id, 95
## 
## Fixed effects:
##               Estimate Std. Error         df t value Pr(>|t|)    
## (Intercept)   6.286355   0.228788 284.644902  27.477   <2e-16 ***
## mri_age_y     0.016356   0.034473 318.479871   0.474    0.635    
## female       -0.003854   0.157311  78.609245  -0.024    0.981    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##           (Intr) mr_g_y
## mri_age_y -0.881       
## female    -0.284 -0.046
anova(lag_ins_age,lag_ins_age_sex)
## Data: lag_dat
## Models:
## lag_ins_age: Ins_conc_molal ~ mri_age_y + (1 | subj_id)
## lag_ins_age_sex: Ins_conc_molal ~ mri_age_y + female + (1 | subj_id)
##                 npar    AIC    BIC  logLik deviance Chisq Df Pr(>Chisq)
## lag_ins_age        4 1013.3 1028.3 -502.64   1005.3                    
## lag_ins_age_sex    5 1015.3 1034.1 -502.64   1005.3 6e-04  1     0.9807
lag_ins_ageXsex<- lmer(Ins_conc_molal ~ mri_age_y + female + mri_age_y:female + (1 | subj_id), data=lag_dat, REML=FALSE)
summary(lag_ins_ageXsex)
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
##   method [lmerModLmerTest]
## Formula: Ins_conc_molal ~ mri_age_y + female + mri_age_y:female + (1 |  
##     subj_id)
##    Data: lag_dat
## 
##      AIC      BIC   logLik deviance df.resid 
##   1017.1   1039.7   -502.6   1005.1      313 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.4248 -0.5418 -0.0121  0.5899  4.6348 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  subj_id  (Intercept) 0.1767   0.4204  
##  Residual             1.2221   1.1055  
## Number of obs: 319, groups:  subj_id, 95
## 
## Fixed effects:
##                    Estimate Std. Error         df t value Pr(>|t|)    
## (Intercept)        6.361751   0.294417 312.731852  21.608   <2e-16 ***
## mri_age_y          0.003458   0.046822 316.620308   0.074    0.941    
## female            -0.171702   0.441468 305.644875  -0.389    0.698    
## mri_age_y:female   0.028161   0.069170 318.668940   0.407    0.684    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) mr_g_y female
## mri_age_y   -0.930              
## female      -0.667  0.620       
## mri_g_y:fml  0.630 -0.677 -0.934
anova(lag_ins_age,lag_ins_ageXsex)
## Data: lag_dat
## Models:
## lag_ins_age: Ins_conc_molal ~ mri_age_y + (1 | subj_id)
## lag_ins_ageXsex: Ins_conc_molal ~ mri_age_y + female + mri_age_y:female + (1 | subj_id)
##                 npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)
## lag_ins_age        4 1013.3 1028.3 -502.64   1005.3                     
## lag_ins_ageXsex    6 1017.1 1039.7 -502.56   1005.1 0.1663  2     0.9202

Glx

lag_glx_age<- lmer(Glu_Gln_conc_molal ~ mri_age_y + (1 | subj_id), data=lag_dat, REML=FALSE)
## boundary (singular) fit: see ?isSingular
isSingular(lag_glx_age) # TRUE: model is almost/near singular (parameters are on the boundary of the feasible parameter space: random effects variance = 0)
## [1] TRUE
summary(lag_glx_age)
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
##   method [lmerModLmerTest]
## Formula: Glu_Gln_conc_molal ~ mri_age_y + (1 | subj_id)
##    Data: lag_dat
## 
##      AIC      BIC   logLik deviance df.resid 
##   1396.1   1411.1   -694.0   1388.1      315 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.6700 -0.6848  0.0127  0.6729  3.5338 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  subj_id  (Intercept) 0.000    0.000   
##  Residual             4.542    2.131   
## Number of obs: 319, groups:  subj_id, 95
## 
## Fixed effects:
##              Estimate Std. Error        df t value Pr(>|t|)    
## (Intercept)  23.23676    0.39326 319.00000  59.088  < 2e-16 ***
## mri_age_y    -0.25961    0.06309 319.00000  -4.115 4.93e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##           (Intr)
## mri_age_y -0.953
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see ?isSingular
lag_glx_age_sex<- lmer(Glu_Gln_conc_molal ~ mri_age_y + female + (1 | subj_id), data=lag_dat, REML=FALSE)
## boundary (singular) fit: see ?isSingular
summary(lag_glx_age_sex)
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
##   method [lmerModLmerTest]
## Formula: Glu_Gln_conc_molal ~ mri_age_y + female + (1 | subj_id)
##    Data: lag_dat
## 
##      AIC      BIC   logLik deviance df.resid 
##   1396.7   1415.5   -693.4   1386.7      314 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.7491 -0.6907  0.0192  0.7055  3.6028 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  subj_id  (Intercept) 0.000    0.000   
##  Residual             4.523    2.127   
## Number of obs: 319, groups:  subj_id, 95
## 
## Fixed effects:
##              Estimate Std. Error        df t value Pr(>|t|)    
## (Intercept)  23.12126    0.40478 319.00000  57.120  < 2e-16 ***
## mri_age_y    -0.26202    0.06299 319.00000  -4.160  4.1e-05 ***
## female        0.27789    0.23880 319.00000   1.164    0.245    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##           (Intr) mr_g_y
## mri_age_y -0.915       
## female    -0.245 -0.033
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see ?isSingular
anova(lag_glx_age,lag_glx_age_sex)
## Data: lag_dat
## Models:
## lag_glx_age: Glu_Gln_conc_molal ~ mri_age_y + (1 | subj_id)
## lag_glx_age_sex: Glu_Gln_conc_molal ~ mri_age_y + female + (1 | subj_id)
##                 npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)
## lag_glx_age        4 1396.1 1411.1 -694.04   1388.1                     
## lag_glx_age_sex    5 1396.7 1415.5 -693.36   1386.7 1.3513  1      0.245
lag_glx_ageXsex<- lmer(Glu_Gln_conc_molal ~ mri_age_y + female + mri_age_y:female + (1 | subj_id), data=lag_dat, REML=FALSE)
## boundary (singular) fit: see ?isSingular
summary(lag_glx_ageXsex)
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
##   method [lmerModLmerTest]
## Formula: Glu_Gln_conc_molal ~ mri_age_y + female + mri_age_y:female +  
##     (1 | subj_id)
##    Data: lag_dat
## 
##      AIC      BIC   logLik deviance df.resid 
##   1397.8   1420.4   -692.9   1385.8      313 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.7059 -0.6942 -0.0013  0.6785  3.6197 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  subj_id  (Intercept) 0.00     0.000   
##  Residual             4.51     2.124   
## Number of obs: 319, groups:  subj_id, 95
## 
## Fixed effects:
##                   Estimate Std. Error        df t value Pr(>|t|)    
## (Intercept)       23.44862    0.52588 319.00000  44.589  < 2e-16 ***
## mri_age_y         -0.31768    0.08502 319.00000  -3.737 0.000221 ***
## female            -0.45362    0.78868 319.00000  -0.575 0.565583    
## mri_age_y:female   0.12295    0.12636 319.00000   0.973 0.331265    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) mr_g_y female
## mri_age_y   -0.951              
## female      -0.667  0.634       
## mri_g_y:fml  0.640 -0.673 -0.953
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see ?isSingular
anova(lag_glx_age,lag_glx_ageXsex)
## Data: lag_dat
## Models:
## lag_glx_age: Glu_Gln_conc_molal ~ mri_age_y + (1 | subj_id)
## lag_glx_ageXsex: Glu_Gln_conc_molal ~ mri_age_y + female + mri_age_y:female + (1 | subj_id)
##                 npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)
## lag_glx_age        4 1396.1 1411.1 -694.04   1388.1                     
## lag_glx_ageXsex    6 1397.8 1420.4 -692.89   1385.8 2.2967  2     0.3172

Plot results LAG

Create scatter plot with regression line based on model predictions NAA

#pull intercept and slope from model for regression line
fixef.lag_naa<-fixef(lag_naa_age)

#plot
plot_lag_naa_age<-ggplot(lag_dat, aes(x=mri_age_y, y=NAA_conc_molal,color=as.factor(female)))
plot_lag_naa_age + 
  geom_point() +
  geom_line(aes(group = as.factor(subj_id))) +
  labs(title = "Left Temporo-Parietal", x = "Age (years)", y = "tNAA (molal)") +
  scale_color_discrete(name = "Sex", labels = c("Male", "Female")) +
  geom_abline(aes(intercept = fixef.lag_naa[[1]], slope= fixef.lag_naa[[2]])) + 
  theme_classic()
## Warning: Removed 3 rows containing missing values (geom_point).
## Warning: Removed 3 row(s) containing missing values (geom_path).

ggsave("/Users/meaghan/OneDrive - University of Calgary/Preschool_data/MRS/results/lag_naa_age.png")
## Saving 7 x 5 in image
## Warning: Removed 3 rows containing missing values (geom_point).

## Warning: Removed 3 row(s) containing missing values (geom_path).

Cho

#pull intercept and slope from model for regression line
fixef.lag_cho<-fixef(lag_cho_age)

#plot
plot_lag_cho_age<-ggplot(lag_dat, aes(x=mri_age_y, y=Cho_GPC_PCh_conc_molal,color=as.factor(female)))
plot_lag_cho_age + 
  geom_point() +
  geom_line(aes(group = as.factor(subj_id))) +
  labs(title = "Left Temporo-Parietal", x = "Age (years)", y = "tCho (molal)") +
  scale_color_discrete(name = "Sex", labels = c("Male", "Female")) +
  geom_abline(aes(intercept = fixef.lag_cho[[1]], slope= fixef.lag_cho[[2]])) + 
  theme_classic()
## Warning: Removed 1 rows containing missing values (geom_point).
## Warning: Removed 1 row(s) containing missing values (geom_path).

ggsave("/Users/meaghan/OneDrive - University of Calgary/Preschool_data/MRS/results/lag_cho_age.png")
## Saving 7 x 5 in image
## Warning: Removed 1 rows containing missing values (geom_point).

## Warning: Removed 1 row(s) containing missing values (geom_path).

Glx

#pull intercept and slope from model for regression line
fixef.lag_glx<-fixef(lag_glx_age)

#plot
plot_lag_glx_age<-ggplot(lag_dat, aes(x=mri_age_y, y=Glu_Gln_conc_molal,color=as.factor(female)))
plot_lag_glx_age + 
  geom_point() +
  geom_line(aes(group = as.factor(subj_id))) +
  labs(title = "Left Temporo-Parietal", x = "Age (years)", y = "Glx (molal)") +
  scale_color_discrete(name = "Sex", labels = c("Male", "Female")) +
  geom_abline(aes(intercept = fixef.lag_glx[[1]], slope= fixef.lag_glx[[2]])) + 
  theme_classic()
## Warning: Removed 2 rows containing missing values (geom_point).
## Warning: Removed 2 row(s) containing missing values (geom_path).

ggsave("/Users/meaghan/OneDrive - University of Calgary/Preschool_data/MRS/results/lag_glx_age.png")
## Saving 7 x 5 in image
## Warning: Removed 2 rows containing missing values (geom_point).

## Warning: Removed 2 row(s) containing missing values (geom_path).

##Correct for multiple comparisons Use FDR method (Benjamini & Hochberg) to correct for mutliple comparisons (5 metabolites x 2 voxels = 10 models)

#list models (voxel/metabolite pair)
models<- c("acg_naa", "acg_cr", "acg_cho", "acg_ins", "acg_glx", "lag_naa", "lag_cr", "lag_cho", "lag_ins", "lag_glx")

#list p-values for each model, matching order of models
p<-c(.0217, .307, .0258, .0494, .23, .00093, .115, .00229, .636, 0)

fdr_table<-data.frame(models, p)

fdr_table$p_fdr<-p.adjust(p, method = "fdr", n= length(p))
fdr_table
##     models       p       p_fdr
## 1  acg_naa 0.02170 0.051600000
## 2   acg_cr 0.30700 0.341111111
## 3  acg_cho 0.02580 0.051600000
## 4  acg_ins 0.04940 0.082333333
## 5  acg_glx 0.23000 0.287500000
## 6  lag_naa 0.00093 0.004650000
## 7   lag_cr 0.11500 0.164285714
## 8  lag_cho 0.00229 0.007633333
## 9  lag_ins 0.63600 0.636000000
## 10 lag_glx 0.00000 0.000000000
#run FDR without LAG-Glx in case it's invalid:
p2<-p<-c(.0217, .307, .0258, .0494, .23, .00093, .115, .00229, .636)
p.adjust(p2, method = "fdr", n= length(p2))
## [1] 0.0580500 0.3453750 0.0580500 0.0889200 0.2957143 0.0083700 0.1725000
## [8] 0.0103050 0.6360000

##Combine plots into 1 figure Simplify plots to legend & region don’t repeat in each panel

library(patchwork)
p1_acg_naa_age<-plot_acg_naa_age + 
  geom_point() +
  geom_line(aes(group = as.factor(subj_id))) +
  labs(title = "Anterior Cingulate", x = "Age (years)", y = "tNAA (molal)") +
  scale_color_discrete(name = "Sex", labels = c("Male", "Female")) +
  geom_abline(aes(intercept = fixef.acg_naa[[1]], slope= fixef.acg_naa[[2]])) + 
  theme_classic() +
  theme(plot.title= element_text(hjust = 0.5))

p4_lag_naa_age<-plot_lag_naa_age + 
  geom_point() +
  geom_line(aes(group = as.factor(subj_id))) +
  labs(title = "Left Temporo-Parietal", x = "Age (years)", y = "tNAA (molal)") +
  scale_color_discrete(guide="legend", name = "Sex", labels = c("Male", "Female")) +
  theme_classic() +
  #theme(legend.position = c(0.7, 0.95), legend.background = element_rect(fill='transparent')) +
  geom_abline(aes(intercept = fixef.lag_naa[[1]], slope= fixef.lag_naa[[2]])) +
  theme(plot.title= element_text(hjust = 0.5))


p2_acg_cho_age<-plot_acg_cho_age + 
  geom_point() +
  geom_line(aes(group = as.factor(subj_id))) +
  labs(x = "Age (years)", y = "tCho (molal)") +
  scale_color_discrete(name = "Sex", labels = c("Male", "Female")) +
  geom_abline(aes(intercept = fixef.acg_cho[[1]], slope= fixef.acg_cho[[2]])) + 
  theme_classic() 
  

p5_lag_cho_age<-plot_lag_cho_age + 
  geom_point() +
  geom_line(aes(group = as.factor(subj_id))) +
  labs(x = "Age (years)", y = "tCho (molal)") +
  scale_color_discrete(name = "Sex", labels = c("Male", "Female")) +
  geom_abline(aes(intercept = fixef.lag_cho[[1]], slope= fixef.lag_cho[[2]])) + 
  theme_classic()

p3_acg_ins_age<-plot_acg_ins_age + 
  geom_point() +
  geom_line(aes(group = as.factor(subj_id))) +
  labs(x = "Age (years)", y = "Ins (molal)") +
  scale_color_discrete(name = "Sex", labels = c("Male", "Female")) +
  geom_abline(aes(intercept = fixef.acg_ins[[1]], slope= fixef.acg_ins[[2]])) + 
  theme_classic()

p6_lag_glx_age<-plot_lag_glx_age + 
  geom_point() +
  geom_line(aes(group = as.factor(subj_id))) +
  labs(x = "Age (years)", y = "Glx (molal)") +
  scale_color_discrete(name = "Sex", labels = c("Male", "Female")) +
  geom_abline(aes(intercept = fixef.lag_glx[[1]], slope= fixef.lag_glx[[2]])) + 
  theme_classic()

results_fig<-p1_acg_naa_age + p4_lag_naa_age + p2_acg_cho_age + p5_lag_cho_age + p3_acg_ins_age + p6_lag_glx_age + plot_layout(nrow=3, guides ='collect')

ggsave("/Users/meaghan/OneDrive - University of Calgary/Preschool_data/MRS/results/results_fig_grid.png", results_fig, width = 2000, units = "px")
## Saving 2000 x 1500 px image
## Warning: Removed 3 rows containing missing values (geom_point).
## Warning: Removed 3 row(s) containing missing values (geom_path).
## Warning: Removed 1 rows containing missing values (geom_point).
## Warning: Removed 1 row(s) containing missing values (geom_path).
## Warning: Removed 1 rows containing missing values (geom_point).
## Warning: Removed 1 row(s) containing missing values (geom_path).
## Warning: Removed 1 rows containing missing values (geom_point).
## Warning: Removed 1 row(s) containing missing values (geom_path).
## Warning: Removed 2 rows containing missing values (geom_point).
## Warning: Removed 2 row(s) containing missing values (geom_path).